{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Packages\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import warnings\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "warnings.filterwarnings('ignore') "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# this installs the XGBoost package\n",
    "# !pip install xgboost      "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# this installs the mlxtend package used for model stacking\n",
    "# !pip install mlxtend      "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.model_selection import GridSearchCV, RandomizedSearchCV\n",
    "from sklearn.metrics import mean_squared_error, r2_score,  mean_absolute_error\n",
    "\n",
    "from sklearn.pipeline import Pipeline\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "\n",
    "from sklearn.linear_model import LinearRegression, LassoCV, RidgeCV, ElasticNetCV\n",
    "from sklearn.tree import DecisionTreeRegressor\n",
    "from sklearn.ensemble import RandomForestRegressor, BaggingRegressor\n",
    "import xgboost as xgb\n",
    "from mlxtend.regressor import StackingCVRegressor\n",
    "from sklearn.tree import export_graphviz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "sns.set_palette([\"#0b559f\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## House Pricing Data\n",
    "\n",
    "The original dataset has been cleaned and processed into a version that is ready for analysis. The dataset has 196 variables after coding the categorical predictors using dummy variables and the creation of other relevant variables. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>1stFlrSF</th>\n",
       "      <th>2ndFlrSF</th>\n",
       "      <th>3SsnPorch</th>\n",
       "      <th>Age</th>\n",
       "      <th>BsmtFinSF1</th>\n",
       "      <th>BsmtFinSF2</th>\n",
       "      <th>BsmtUnfSF</th>\n",
       "      <th>EnclosedPorch</th>\n",
       "      <th>GarageArea</th>\n",
       "      <th>LotArea</th>\n",
       "      <th>...</th>\n",
       "      <th>RoofMatl_Other</th>\n",
       "      <th>RoofStyle_Hip</th>\n",
       "      <th>RoofStyle_Other</th>\n",
       "      <th>ScreenPorchZero</th>\n",
       "      <th>WoodDeckSFZero</th>\n",
       "      <th>YrSold_2007</th>\n",
       "      <th>YrSold_2008</th>\n",
       "      <th>YrSold_2009</th>\n",
       "      <th>YrSold_2010</th>\n",
       "      <th>SalePrice</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1656</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>50</td>\n",
       "      <td>639.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>441.0</td>\n",
       "      <td>0</td>\n",
       "      <td>528.0</td>\n",
       "      <td>31770</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>215000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>896</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>49</td>\n",
       "      <td>468.0</td>\n",
       "      <td>144.0</td>\n",
       "      <td>270.0</td>\n",
       "      <td>0</td>\n",
       "      <td>730.0</td>\n",
       "      <td>11622</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>105000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1329</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>52</td>\n",
       "      <td>923.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>406.0</td>\n",
       "      <td>0</td>\n",
       "      <td>312.0</td>\n",
       "      <td>14267</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>172000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>2110</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>42</td>\n",
       "      <td>1065.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1045.0</td>\n",
       "      <td>0</td>\n",
       "      <td>522.0</td>\n",
       "      <td>11160</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>244000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>928</td>\n",
       "      <td>701</td>\n",
       "      <td>0</td>\n",
       "      <td>13</td>\n",
       "      <td>791.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>137.0</td>\n",
       "      <td>0</td>\n",
       "      <td>482.0</td>\n",
       "      <td>13830</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>189900</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 196 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "   1stFlrSF  2ndFlrSF  3SsnPorch  Age  BsmtFinSF1  BsmtFinSF2  BsmtUnfSF  \\\n",
       "0      1656         0          0   50       639.0         0.0      441.0   \n",
       "1       896         0          0   49       468.0       144.0      270.0   \n",
       "2      1329         0          0   52       923.0         0.0      406.0   \n",
       "3      2110         0          0   42      1065.0         0.0     1045.0   \n",
       "4       928       701          0   13       791.0         0.0      137.0   \n",
       "\n",
       "   EnclosedPorch  GarageArea  LotArea  ...  RoofMatl_Other  RoofStyle_Hip  \\\n",
       "0              0       528.0    31770  ...               0              1   \n",
       "1              0       730.0    11622  ...               0              0   \n",
       "2              0       312.0    14267  ...               0              1   \n",
       "3              0       522.0    11160  ...               0              1   \n",
       "4              0       482.0    13830  ...               0              0   \n",
       "\n",
       "   RoofStyle_Other  ScreenPorchZero  WoodDeckSFZero  YrSold_2007  YrSold_2008  \\\n",
       "0                0                1               0            0            0   \n",
       "1                0                0               0            0            0   \n",
       "2                0                1               0            0            0   \n",
       "3                0                1               1            0            0   \n",
       "4                0                1               0            0            0   \n",
       "\n",
       "   YrSold_2009  YrSold_2010  SalePrice  \n",
       "0            0            1     215000  \n",
       "1            0            1     105000  \n",
       "2            0            1     172000  \n",
       "3            0            1     244000  \n",
       "4            0            1     189900  \n",
       "\n",
       "[5 rows x 196 columns]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data=pd.read_csv('Datasets/AmesHousing-Processed.csv')\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As usual, we the split the data into training (70%) and test sets (30%)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "response='SalePrice'\n",
    "predictors=list(data.columns.values[:-1])\n",
    "\n",
    "# Randomly split indexes\n",
    "index_train, index_test  = train_test_split(np.array(data.index), train_size=0.7, random_state=5)\n",
    "\n",
    "# Write training and test sets \n",
    "train = data.loc[index_train,:].copy()\n",
    "test =  data.loc[index_test,:].copy()\n",
    "\n",
    "# Write training and test response vectors\n",
    "y_train = np.log(train[response])\n",
    "y_test = np.log(test[response])\n",
    "\n",
    "# Write training and test design matrices\n",
    "X_train = train[predictors].copy()\n",
    "X_test = test[predictors].copy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that our Y is the logarithm of the Sale Price."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we implement a number of familiar methods."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear Regression "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LinearRegression(copy_X=True, fit_intercept=True, n_jobs=None, normalize=False)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ols = LinearRegression()\n",
    "ols.fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Regularised Linear Models\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Lasso"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1.3 s, sys: 13.3 ms, total: 1.31 s\n",
      "Wall time: 330 ms\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Pipeline(memory=None,\n",
       "         steps=[('scaler',\n",
       "                 StandardScaler(copy=True, with_mean=True, with_std=True)),\n",
       "                ('estimator',\n",
       "                 LassoCV(alphas=None, copy_X=True, cv=5, eps=0.001,\n",
       "                         fit_intercept=True, max_iter=1000, n_alphas=100,\n",
       "                         n_jobs=None, normalize=False, positive=False,\n",
       "                         precompute='auto', random_state=None,\n",
       "                         selection='cyclic', tol=0.0001, verbose=False))],\n",
       "         verbose=False)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "lasso = Pipeline((\n",
    "    ('scaler', StandardScaler()),\n",
    "    ('estimator', LassoCV(cv=5)),\n",
    "))\n",
    "\n",
    "lasso.fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ridge Regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 10.2 s, sys: 61.7 ms, total: 10.3 s\n",
      "Wall time: 2.58 s\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Pipeline(memory=None,\n",
       "         steps=[('scaler',\n",
       "                 StandardScaler(copy=True, with_mean=True, with_std=True)),\n",
       "                ('estimator',\n",
       "                 RidgeCV(alphas=array([3.05175781e-05, 3.50554918e-05, 4.02681858e-05, 4.62559987e-05,\n",
       "       5.31341897e-05, 6.10351562e-05, 7.01109836e-05, 8.05363715e-05,\n",
       "       9.25119975e-05, 1.06268379e-04, 1.22070312e-04, 1.40221967e-04,\n",
       "       1.61072743e-04, 1.85023995e-04, 2.12536759e...\n",
       "       2.70235220e+03, 3.10418753e+03, 3.56577511e+03, 4.09600000e+03,\n",
       "       4.70506846e+03, 5.40470440e+03, 6.20837506e+03, 7.13155021e+03,\n",
       "       8.19200000e+03, 9.41013692e+03, 1.08094088e+04, 1.24167501e+04,\n",
       "       1.42631004e+04, 1.63840000e+04, 1.88202738e+04, 2.16188176e+04,\n",
       "       2.48335002e+04, 2.85262009e+04, 3.27680000e+04]),\n",
       "                         cv=5, fit_intercept=True, gcv_mode=None,\n",
       "                         normalize=False, scoring=None,\n",
       "                         store_cv_values=False))],\n",
       "         verbose=False)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "alphas = list(np.logspace(-15, 15, 151, base=2))\n",
    "\n",
    "ridge = Pipeline((\n",
    "    ('scaler', StandardScaler()),\n",
    "    ('estimator', RidgeCV(alphas=alphas, cv=5)),\n",
    "))\n",
    "\n",
    "ridge.fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Elastic Net"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 12.6 s, sys: 57.7 ms, total: 12.6 s\n",
      "Wall time: 3.17 s\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Pipeline(memory=None,\n",
       "         steps=[('scaler',\n",
       "                 StandardScaler(copy=True, with_mean=True, with_std=True)),\n",
       "                ('estimator',\n",
       "                 ElasticNetCV(alphas=None, copy_X=True, cv=5, eps=0.001,\n",
       "                              fit_intercept=True,\n",
       "                              l1_ratio=[0.01, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7,\n",
       "                                        0.8, 0.9, 0.99],\n",
       "                              max_iter=1000, n_alphas=100, n_jobs=None,\n",
       "                              normalize=False, positive=False,\n",
       "                              precompute='auto', random_state=None,\n",
       "                              selection='cyclic', tol=0.0001, verbose=0))],\n",
       "         verbose=False)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "enet = Pipeline((\n",
    "    ('scaler', StandardScaler()),\n",
    "    ('estimator', ElasticNetCV(l1_ratio=[0.01,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9, 0.99], cv=5)),\n",
    "))\n",
    "\n",
    "enet.fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Regression Tree"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Best parameters: {'min_samples_leaf': 5, 'max_depth': 17}\n",
      "CPU times: user 2.27 s, sys: 28.5 ms, total: 2.3 s\n",
      "Wall time: 1.76 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "model = DecisionTreeRegressor(min_samples_leaf=5)\n",
    "\n",
    "tuning_parameters = {\n",
    "    'min_samples_leaf': [1,5,10,20],\n",
    "    'max_depth': np.arange(1,30),\n",
    "}\n",
    "\n",
    "tree = RandomizedSearchCV(model, tuning_parameters, n_iter=20, cv=5, return_train_score=False)\n",
    "tree.fit(X_train, y_train)\n",
    "\n",
    "print('Best parameters:', tree.best_params_)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we visualise the output from the regression tree. Please note, we are constructing an artificially simple model for the sake of interpretability. In fact, one of the benefits of decision trees is their ability to construct large and often complicated decision paths. Looking at the difference between the maximum depth (3 vs 30) and the minimum samples in the terminal nodes (150 vs 1) highlights the corresponding reduction in the size of the tree. Although the visualised tree model is not as complex, this visualisation is still insightful."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "digraph Tree {\n",
      "node [shape=box] ;\n",
      "0 [label=\"X[3] <= 23.5\\nmse = 0.141\\nsamples = 1688\\nvalue = 12.008\"] ;\n",
      "1 [label=\"X[25] <= 2.5\\nmse = 0.079\\nsamples = 614\\nvalue = 12.311\"] ;\n",
      "0 -> 1 [labeldistance=2.5, labelangle=45, headlabel=\"True\"] ;\n",
      "2 [label=\"X[0] <= 1270.5\\nmse = 0.042\\nsamples = 457\\nvalue = 12.205\"] ;\n",
      "1 -> 2 ;\n",
      "3 [label=\"mse = 0.031\\nsamples = 302\\nvalue = 12.136\"] ;\n",
      "2 -> 3 ;\n",
      "4 [label=\"mse = 0.038\\nsamples = 155\\nvalue = 12.339\"] ;\n",
      "2 -> 4 ;\n",
      "5 [label=\"mse = 0.059\\nsamples = 157\\nvalue = 12.62\"] ;\n",
      "1 -> 5 ;\n",
      "6 [label=\"X[86] <= 0.5\\nmse = 0.094\\nsamples = 1074\\nvalue = 11.835\"] ;\n",
      "0 -> 6 [labeldistance=2.5, labelangle=-45, headlabel=\"False\"] ;\n",
      "7 [label=\"X[9] <= 10220.0\\nmse = 0.078\\nsamples = 464\\nvalue = 12.008\"] ;\n",
      "6 -> 7 ;\n",
      "8 [label=\"mse = 0.048\\nsamples = 233\\nvalue = 11.871\"] ;\n",
      "7 -> 8 ;\n",
      "9 [label=\"mse = 0.07\\nsamples = 231\\nvalue = 12.147\"] ;\n",
      "7 -> 9 ;\n",
      "10 [label=\"X[18] <= 682.0\\nmse = 0.066\\nsamples = 610\\nvalue = 11.702\"] ;\n",
      "6 -> 10 ;\n",
      "11 [label=\"mse = 0.085\\nsamples = 158\\nvalue = 11.495\"] ;\n",
      "10 -> 11 ;\n",
      "12 [label=\"mse = 0.039\\nsamples = 452\\nvalue = 11.775\"] ;\n",
      "10 -> 12 ;\n",
      "}\n"
     ]
    }
   ],
   "source": [
    "# Please note how much simpler the model we plot is when compared to the actual model. \n",
    "# For example, look at the maximum depth \n",
    "from sklearn.tree import export_graphviz\n",
    "model_plot = DecisionTreeRegressor(max_depth=3, min_samples_leaf=150)\n",
    "tuning_parameters_plot = {\n",
    "    'min_samples_leaf': [150],\n",
    "    'max_depth': [1,2,3],\n",
    "}\n",
    "tree_plot = RandomizedSearchCV(model_plot, tuning_parameters_plot, n_iter=20, cv=5, return_train_score=False)\n",
    "tree_plot.fit(X_train, y_train)\n",
    "best_estimator_plot = tree_plot.best_estimator_\n",
    "dot_data = export_graphviz(best_estimator_plot)\n",
    "print(dot_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 2.40.1 (20161225.0304)\n",
       " -->\n",
       "<!-- Title: Tree Pages: 1 -->\n",
       "<svg width=\"803pt\" height=\"313pt\"\n",
       " viewBox=\"0.00 0.00 803.00 313.00\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 309)\">\n",
       "<title>Tree</title>\n",
       "<polygon fill=\"#ffffff\" stroke=\"transparent\" points=\"-4,4 -4,-309 799,-309 799,4 -4,4\"/>\n",
       "<!-- 0 -->\n",
       "<g id=\"node1\" class=\"node\">\n",
       "<title>0</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M422.5,-305C422.5,-305 330.5,-305 330.5,-305 324.5,-305 318.5,-299 318.5,-293 318.5,-293 318.5,-264 318.5,-264 318.5,-258 324.5,-252 330.5,-252 330.5,-252 422.5,-252 422.5,-252 428.5,-252 434.5,-258 434.5,-264 434.5,-264 434.5,-293 434.5,-293 434.5,-299 428.5,-305 422.5,-305\"/>\n",
       "<text text-anchor=\"middle\" x=\"376.5\" y=\"-289.8\" font-family=\"Helvetica,sans-Serif\" font-size=\"14.00\" fill=\"#000000\">Age &lt;= 23.5</text>\n",
       "<text text-anchor=\"middle\" x=\"376.5\" y=\"-274.8\" font-family=\"Helvetica,sans-Serif\" font-size=\"14.00\" fill=\"#000000\">samples = 1688</text>\n",
       "<text text-anchor=\"middle\" x=\"376.5\" y=\"-259.8\" font-family=\"Helvetica,sans-Serif\" font-size=\"14.00\" fill=\"#000000\">value = 12.008</text>\n",
       "</g>\n",
       "<!-- 1 -->\n",
       "<g id=\"node2\" class=\"node\">\n",
       "<title>1</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M341.5,-216C341.5,-216 229.5,-216 229.5,-216 223.5,-216 217.5,-210 217.5,-204 217.5,-204 217.5,-175 217.5,-175 217.5,-169 223.5,-163 229.5,-163 229.5,-163 341.5,-163 341.5,-163 347.5,-163 353.5,-169 353.5,-175 353.5,-175 353.5,-204 353.5,-204 353.5,-210 347.5,-216 341.5,-216\"/>\n",
       "<text text-anchor=\"middle\" x=\"285.5\" y=\"-200.8\" font-family=\"Helvetica,sans-Serif\" font-size=\"14.00\" fill=\"#000000\">GarageCars &lt;= 2.5</text>\n",
       "<text text-anchor=\"middle\" x=\"285.5\" y=\"-185.8\" font-family=\"Helvetica,sans-Serif\" font-size=\"14.00\" fill=\"#000000\">samples = 614</text>\n",
       "<text text-anchor=\"middle\" x=\"285.5\" y=\"-170.8\" font-family=\"Helvetica,sans-Serif\" font-size=\"14.00\" fill=\"#000000\">value = 12.311</text>\n",
       "</g>\n",
       "<!-- 0&#45;&gt;1 -->\n",
       "<g id=\"edge1\" class=\"edge\">\n",
       "<title>0&#45;&gt;1</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M349.197,-251.7971C340.0189,-242.8207 329.6581,-232.6876 320.0006,-223.2423\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"322.3802,-220.674 312.7837,-216.1841 317.4857,-225.6784 322.3802,-220.674\"/>\n",
       "<text text-anchor=\"middle\" x=\"313.0615\" y=\"-237.4825\" font-family=\"Helvetica,sans-Serif\" font-size=\"14.00\" fill=\"#000000\">True</text>\n",
       "</g>\n",
       "<!-- 6 -->\n",
       "<g id=\"node7\" class=\"node\">\n",
       "<title>6</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M551.5,-216C551.5,-216 383.5,-216 383.5,-216 377.5,-216 371.5,-210 371.5,-204 371.5,-204 371.5,-175 371.5,-175 371.5,-169 377.5,-163 383.5,-163 383.5,-163 551.5,-163 551.5,-163 557.5,-163 563.5,-169 563.5,-175 563.5,-175 563.5,-204 563.5,-204 563.5,-210 557.5,-216 551.5,-216\"/>\n",
       "<text text-anchor=\"middle\" x=\"467.5\" y=\"-200.8\" font-family=\"Helvetica,sans-Serif\" font-size=\"14.00\" fill=\"#000000\">FireplaceQu_Missing &lt;= 0.5</text>\n",
       "<text text-anchor=\"middle\" x=\"467.5\" y=\"-185.8\" font-family=\"Helvetica,sans-Serif\" font-size=\"14.00\" fill=\"#000000\">samples = 1074</text>\n",
       "<text text-anchor=\"middle\" x=\"467.5\" y=\"-170.8\" font-family=\"Helvetica,sans-Serif\" font-size=\"14.00\" fill=\"#000000\">value = 11.835</text>\n",
       "</g>\n",
       "<!-- 0&#45;&gt;6 -->\n",
       "<g id=\"edge6\" class=\"edge\">\n",
       "<title>0&#45;&gt;6</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M403.803,-251.7971C412.9811,-242.8207 423.3419,-232.6876 432.9994,-223.2423\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"435.5143,-225.6784 440.2163,-216.1841 430.6198,-220.674 435.5143,-225.6784\"/>\n",
       "<text text-anchor=\"middle\" x=\"439.9385\" y=\"-237.4825\" font-family=\"Helvetica,sans-Serif\" font-size=\"14.00\" fill=\"#000000\">False</text>\n",
       "</g>\n",
       "<!-- 2 -->\n",
       "<g id=\"node3\" class=\"node\">\n",
       "<title>2</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M211.5,-127C211.5,-127 99.5,-127 99.5,-127 93.5,-127 87.5,-121 87.5,-115 87.5,-115 87.5,-86 87.5,-86 87.5,-80 93.5,-74 99.5,-74 99.5,-74 211.5,-74 211.5,-74 217.5,-74 223.5,-80 223.5,-86 223.5,-86 223.5,-115 223.5,-115 223.5,-121 217.5,-127 211.5,-127\"/>\n",
       "<text text-anchor=\"middle\" x=\"155.5\" y=\"-111.8\" font-family=\"Helvetica,sans-Serif\" font-size=\"14.00\" fill=\"#000000\">1stFlrSF &lt;= 1270.5</text>\n",
       "<text text-anchor=\"middle\" x=\"155.5\" y=\"-96.8\" font-family=\"Helvetica,sans-Serif\" font-size=\"14.00\" fill=\"#000000\">samples = 457</text>\n",
       "<text text-anchor=\"middle\" x=\"155.5\" y=\"-81.8\" font-family=\"Helvetica,sans-Serif\" font-size=\"14.00\" fill=\"#000000\">value = 12.205</text>\n",
       "</g>\n",
       "<!-- 1&#45;&gt;2 -->\n",
       "<g id=\"edge2\" class=\"edge\">\n",
       "<title>1&#45;&gt;2</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M246.4957,-162.7971C232.7414,-153.3807 217.1278,-142.6914 202.7652,-132.8585\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"204.7055,-129.9452 194.4768,-127.1841 200.7511,-135.7213 204.7055,-129.9452\"/>\n",
       "</g>\n",
       "<!-- 5 -->\n",
       "<g id=\"node6\" class=\"node\">\n",
       "<title>5</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M337.5,-119.5C337.5,-119.5 253.5,-119.5 253.5,-119.5 247.5,-119.5 241.5,-113.5 241.5,-107.5 241.5,-107.5 241.5,-93.5 241.5,-93.5 241.5,-87.5 247.5,-81.5 253.5,-81.5 253.5,-81.5 337.5,-81.5 337.5,-81.5 343.5,-81.5 349.5,-87.5 349.5,-93.5 349.5,-93.5 349.5,-107.5 349.5,-107.5 349.5,-113.5 343.5,-119.5 337.5,-119.5\"/>\n",
       "<text text-anchor=\"middle\" x=\"295.5\" y=\"-104.3\" font-family=\"Helvetica,sans-Serif\" font-size=\"14.00\" fill=\"#000000\">samples = 157</text>\n",
       "<text text-anchor=\"middle\" x=\"295.5\" y=\"-89.3\" font-family=\"Helvetica,sans-Serif\" font-size=\"14.00\" fill=\"#000000\">value = 12.62</text>\n",
       "</g>\n",
       "<!-- 1&#45;&gt;5 -->\n",
       "<g id=\"edge5\" class=\"edge\">\n",
       "<title>1&#45;&gt;5</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M288.5003,-162.7971C289.674,-152.3518 291.0236,-140.3403 292.2226,-129.6689\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"295.7208,-129.8801 293.3594,-119.5518 288.7646,-129.0984 295.7208,-129.8801\"/>\n",
       "</g>\n",
       "<!-- 3 -->\n",
       "<g id=\"node4\" class=\"node\">\n",
       "<title>3</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M97,-38C97,-38 12,-38 12,-38 6,-38 0,-32 0,-26 0,-26 0,-12 0,-12 0,-6 6,0 12,0 12,0 97,0 97,0 103,0 109,-6 109,-12 109,-12 109,-26 109,-26 109,-32 103,-38 97,-38\"/>\n",
       "<text text-anchor=\"middle\" x=\"54.5\" y=\"-22.8\" font-family=\"Helvetica,sans-Serif\" font-size=\"14.00\" fill=\"#000000\">samples = 302</text>\n",
       "<text text-anchor=\"middle\" x=\"54.5\" y=\"-7.8\" font-family=\"Helvetica,sans-Serif\" font-size=\"14.00\" fill=\"#000000\">value = 12.136</text>\n",
       "</g>\n",
       "<!-- 2&#45;&gt;3 -->\n",
       "<g id=\"edge3\" class=\"edge\">\n",
       "<title>2&#45;&gt;3</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M122.4288,-73.8139C110.8552,-64.4747 97.8895,-54.0123 86.4124,-44.7511\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"88.3399,-41.8091 78.3596,-38.2531 83.944,-47.2567 88.3399,-41.8091\"/>\n",
       "</g>\n",
       "<!-- 4 -->\n",
       "<g id=\"node5\" class=\"node\">\n",
       "<title>4</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M224,-38C224,-38 139,-38 139,-38 133,-38 127,-32 127,-26 127,-26 127,-12 127,-12 127,-6 133,0 139,0 139,0 224,0 224,0 230,0 236,-6 236,-12 236,-12 236,-26 236,-26 236,-32 230,-38 224,-38\"/>\n",
       "<text text-anchor=\"middle\" x=\"181.5\" y=\"-22.8\" font-family=\"Helvetica,sans-Serif\" font-size=\"14.00\" fill=\"#000000\">samples = 155</text>\n",
       "<text text-anchor=\"middle\" x=\"181.5\" y=\"-7.8\" font-family=\"Helvetica,sans-Serif\" font-size=\"14.00\" fill=\"#000000\">value = 12.339</text>\n",
       "</g>\n",
       "<!-- 2&#45;&gt;4 -->\n",
       "<g id=\"edge4\" class=\"edge\">\n",
       "<title>2&#45;&gt;4</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M164.0134,-73.8139C166.6456,-65.5628 169.5576,-56.4349 172.2378,-48.0333\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"175.653,-48.8438 175.3579,-38.2531 168.9841,-46.7163 175.653,-48.8438\"/>\n",
       "</g>\n",
       "<!-- 7 -->\n",
       "<g id=\"node8\" class=\"node\">\n",
       "<title>7</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M515.5,-127C515.5,-127 399.5,-127 399.5,-127 393.5,-127 387.5,-121 387.5,-115 387.5,-115 387.5,-86 387.5,-86 387.5,-80 393.5,-74 399.5,-74 399.5,-74 515.5,-74 515.5,-74 521.5,-74 527.5,-80 527.5,-86 527.5,-86 527.5,-115 527.5,-115 527.5,-121 521.5,-127 515.5,-127\"/>\n",
       "<text text-anchor=\"middle\" x=\"457.5\" y=\"-111.8\" font-family=\"Helvetica,sans-Serif\" font-size=\"14.00\" fill=\"#000000\">LotArea &lt;= 10220.0</text>\n",
       "<text text-anchor=\"middle\" x=\"457.5\" y=\"-96.8\" font-family=\"Helvetica,sans-Serif\" font-size=\"14.00\" fill=\"#000000\">samples = 464</text>\n",
       "<text text-anchor=\"middle\" x=\"457.5\" y=\"-81.8\" font-family=\"Helvetica,sans-Serif\" font-size=\"14.00\" fill=\"#000000\">value = 12.008</text>\n",
       "</g>\n",
       "<!-- 6&#45;&gt;7 -->\n",
       "<g id=\"edge7\" class=\"edge\">\n",
       "<title>6&#45;&gt;7</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M464.4997,-162.7971C463.5999,-154.7887 462.5966,-145.8597 461.6371,-137.3198\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"465.093,-136.7307 460.4982,-127.1841 458.1368,-137.5124 465.093,-136.7307\"/>\n",
       "</g>\n",
       "<!-- 10 -->\n",
       "<g id=\"node11\" class=\"node\">\n",
       "<title>10</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M691,-127C691,-127 558,-127 558,-127 552,-127 546,-121 546,-115 546,-115 546,-86 546,-86 546,-80 552,-74 558,-74 558,-74 691,-74 691,-74 697,-74 703,-80 703,-86 703,-86 703,-115 703,-115 703,-121 697,-127 691,-127\"/>\n",
       "<text text-anchor=\"middle\" x=\"624.5\" y=\"-111.8\" font-family=\"Helvetica,sans-Serif\" font-size=\"14.00\" fill=\"#000000\">TotalBsmtSF &lt;= 682.0</text>\n",
       "<text text-anchor=\"middle\" x=\"624.5\" y=\"-96.8\" font-family=\"Helvetica,sans-Serif\" font-size=\"14.00\" fill=\"#000000\">samples = 610</text>\n",
       "<text text-anchor=\"middle\" x=\"624.5\" y=\"-81.8\" font-family=\"Helvetica,sans-Serif\" font-size=\"14.00\" fill=\"#000000\">value = 11.702</text>\n",
       "</g>\n",
       "<!-- 6&#45;&gt;10 -->\n",
       "<g id=\"edge10\" class=\"edge\">\n",
       "<title>6&#45;&gt;10</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M514.6052,-162.7971C531.6035,-153.1611 550.9533,-142.1921 568.6293,-132.1719\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"570.7254,-135.007 577.6988,-127.0306 567.2733,-128.9174 570.7254,-135.007\"/>\n",
       "</g>\n",
       "<!-- 8 -->\n",
       "<g id=\"node9\" class=\"node\">\n",
       "<title>8</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M385,-38C385,-38 300,-38 300,-38 294,-38 288,-32 288,-26 288,-26 288,-12 288,-12 288,-6 294,0 300,0 300,0 385,0 385,0 391,0 397,-6 397,-12 397,-12 397,-26 397,-26 397,-32 391,-38 385,-38\"/>\n",
       "<text text-anchor=\"middle\" x=\"342.5\" y=\"-22.8\" font-family=\"Helvetica,sans-Serif\" font-size=\"14.00\" fill=\"#000000\">samples = 233</text>\n",
       "<text text-anchor=\"middle\" x=\"342.5\" y=\"-7.8\" font-family=\"Helvetica,sans-Serif\" font-size=\"14.00\" fill=\"#000000\">value = 11.871</text>\n",
       "</g>\n",
       "<!-- 7&#45;&gt;8 -->\n",
       "<g id=\"edge8\" class=\"edge\">\n",
       "<title>7&#45;&gt;8</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M419.8447,-73.8139C406.4109,-64.2934 391.33,-53.6056 378.0766,-44.213\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"379.8495,-41.1796 369.6669,-38.2531 375.802,-46.8909 379.8495,-41.1796\"/>\n",
       "</g>\n",
       "<!-- 9 -->\n",
       "<g id=\"node10\" class=\"node\">\n",
       "<title>9</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M512,-38C512,-38 427,-38 427,-38 421,-38 415,-32 415,-26 415,-26 415,-12 415,-12 415,-6 421,0 427,0 427,0 512,0 512,0 518,0 524,-6 524,-12 524,-12 524,-26 524,-26 524,-32 518,-38 512,-38\"/>\n",
       "<text text-anchor=\"middle\" x=\"469.5\" y=\"-22.8\" font-family=\"Helvetica,sans-Serif\" font-size=\"14.00\" fill=\"#000000\">samples = 231</text>\n",
       "<text text-anchor=\"middle\" x=\"469.5\" y=\"-7.8\" font-family=\"Helvetica,sans-Serif\" font-size=\"14.00\" fill=\"#000000\">value = 12.147</text>\n",
       "</g>\n",
       "<!-- 7&#45;&gt;9 -->\n",
       "<g id=\"edge9\" class=\"edge\">\n",
       "<title>7&#45;&gt;9</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M461.4292,-73.8139C462.6308,-65.6535 463.9586,-56.6354 465.1843,-48.3105\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"468.6711,-48.6563 466.6652,-38.2531 461.7458,-47.6365 468.6711,-48.6563\"/>\n",
       "</g>\n",
       "<!-- 11 -->\n",
       "<g id=\"node12\" class=\"node\">\n",
       "<title>11</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M656,-38C656,-38 571,-38 571,-38 565,-38 559,-32 559,-26 559,-26 559,-12 559,-12 559,-6 565,0 571,0 571,0 656,0 656,0 662,0 668,-6 668,-12 668,-12 668,-26 668,-26 668,-32 662,-38 656,-38\"/>\n",
       "<text text-anchor=\"middle\" x=\"613.5\" y=\"-22.8\" font-family=\"Helvetica,sans-Serif\" font-size=\"14.00\" fill=\"#000000\">samples = 158</text>\n",
       "<text text-anchor=\"middle\" x=\"613.5\" y=\"-7.8\" font-family=\"Helvetica,sans-Serif\" font-size=\"14.00\" fill=\"#000000\">value = 11.495</text>\n",
       "</g>\n",
       "<!-- 10&#45;&gt;11 -->\n",
       "<g id=\"edge11\" class=\"edge\">\n",
       "<title>10&#45;&gt;11</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M620.8982,-73.8139C619.7968,-65.6535 618.5796,-56.6354 617.456,-48.3105\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"620.9048,-47.695 616.0986,-38.2531 613.9677,-48.6314 620.9048,-47.695\"/>\n",
       "</g>\n",
       "<!-- 12 -->\n",
       "<g id=\"node13\" class=\"node\">\n",
       "<title>12</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M783,-38C783,-38 698,-38 698,-38 692,-38 686,-32 686,-26 686,-26 686,-12 686,-12 686,-6 692,0 698,0 698,0 783,0 783,0 789,0 795,-6 795,-12 795,-12 795,-26 795,-26 795,-32 789,-38 783,-38\"/>\n",
       "<text text-anchor=\"middle\" x=\"740.5\" y=\"-22.8\" font-family=\"Helvetica,sans-Serif\" font-size=\"14.00\" fill=\"#000000\">samples = 452</text>\n",
       "<text text-anchor=\"middle\" x=\"740.5\" y=\"-7.8\" font-family=\"Helvetica,sans-Serif\" font-size=\"14.00\" fill=\"#000000\">value = 11.775</text>\n",
       "</g>\n",
       "<!-- 10&#45;&gt;12 -->\n",
       "<g id=\"edge12\" class=\"edge\">\n",
       "<title>10&#45;&gt;12</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M662.4827,-73.8139C676.0333,-64.2934 691.2454,-53.6056 704.614,-44.213\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"706.9266,-46.8658 713.0969,-38.2531 702.9024,-41.1381 706.9266,-46.8658\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.files.Source at 0x7fd378eea050>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import graphviz\n",
    "from sklearn.tree import export_graphviz\n",
    "\n",
    "dot_data = export_graphviz(best_estimator_plot, out_file=None, feature_names=predictors, impurity=False, rounded=True) \n",
    "graph = graphviz.Source(dot_data)\n",
    "graph.render('tree01') # saves tree to a file\n",
    "graph"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bagging\n",
    "\n",
    "We can use the following syntax to implement bagging for the regression trees (however, note that bagging is also a special case of random forests)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 11.4 s, sys: 57.1 ms, total: 11.4 s\n",
      "Wall time: 11.4 s\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "BaggingRegressor(base_estimator=None, bootstrap=True, bootstrap_features=False,\n",
       "                 max_features=1.0, max_samples=1.0, n_estimators=500,\n",
       "                 n_jobs=None, oob_score=False, random_state=1, verbose=0,\n",
       "                 warm_start=False)"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "bag = BaggingRegressor(n_estimators=500, random_state=1)\n",
    "bag.fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Random Forest Regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Best parameters found by randomised search: {'min_samples_leaf': 1, 'max_features': 56} \n",
      "\n",
      "CPU times: user 1.15 s, sys: 79.7 ms, total: 1.23 s\n",
      "Wall time: 18.7 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "model = RandomForestRegressor(n_estimators=100)\n",
    "\n",
    "tuning_parameters = {\n",
    "    'min_samples_leaf': [1,5, 10, 20, 50],\n",
    "    'max_features': np.arange(1, X_train.shape[1], 5),\n",
    "}\n",
    "\n",
    "rf_search = RandomizedSearchCV(model, tuning_parameters, cv = 5, n_iter= 16, return_train_score=False, n_jobs=4,\n",
    "                              random_state = 20)\n",
    "rf_search.fit(X_train, y_train)\n",
    "\n",
    "rf = rf_search.best_estimator_\n",
    "\n",
    "print('Best parameters found by randomised search:', rf_search.best_params_, '\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "RandomForestRegressor(bootstrap=True, ccp_alpha=0.0, criterion='mse',\n",
       "                      max_depth=None, max_features=56, max_leaf_nodes=None,\n",
       "                      max_samples=None, min_impurity_decrease=0.0,\n",
       "                      min_impurity_split=None, min_samples_leaf=1,\n",
       "                      min_samples_split=2, min_weight_fraction_leaf=0.0,\n",
       "                      n_estimators=500, n_jobs=None, oob_score=False,\n",
       "                      random_state=None, verbose=0, warm_start=False)"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rf.n_estimators = 500\n",
    "rf.fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Boosting\n",
    "\n",
    "We start with the <TT>Sciki-Learn</TT> implementation of boosting available in the [<TT>GradientBoostingRegressor</TT>](http://scikit-learn.org/stable/modules/generated/sklearn.ensemble.GradientBoostingRegressor.html) class. Recall that boosting has three crucial tuning parameters:\n",
    "\n",
    "<li style=\"margin-top:15px; margin-bottom: 10px\"> The learning rate.</li> \n",
    "\n",
    "<li style=\"margin-top:10px; margin-bottom: 10px\"> The number of trees.</li> \n",
    "\n",
    "<li style=\"margin-top:10px; margin-bottom: 10px\"> The size of each tree.</li> \n",
    "\n",
    "In addition, we may want to use stochastic gradient boosting by fitting each tree based on a subsample of the training data. \n",
    "\n",
    "The basic syntax for fitting a gradient boosting regressor with <TT>Scikit-Learn</TT> is as follows. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "GradientBoostingRegressor(alpha=0.9, ccp_alpha=0.0, criterion='friedman_mse',\n",
       "                          init=None, learning_rate=0.05, loss='ls', max_depth=4,\n",
       "                          max_features=None, max_leaf_nodes=None,\n",
       "                          min_impurity_decrease=0.0, min_impurity_split=None,\n",
       "                          min_samples_leaf=1, min_samples_split=2,\n",
       "                          min_weight_fraction_leaf=0.0, n_estimators=750,\n",
       "                          n_iter_no_change=None, presort='deprecated',\n",
       "                          random_state=None, subsample=1.0, tol=0.0001,\n",
       "                          validation_fraction=0.1, verbose=0, warm_start=False)"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn.ensemble import GradientBoostingRegressor\n",
    "\n",
    "gb = GradientBoostingRegressor(learning_rate= 0.05, max_depth = 4, n_estimators= 750, subsample = 1.0)\n",
    "gb.fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use a randomised search to tune the model. Note that it useful to keep track of the running time, as the presence of multiple tuning parameters can make this process slow. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Best parameters found by randomised search: {'subsample': 1.0, 'n_estimators': 500, 'max_depth': 2, 'learning_rate': 0.01} \n",
      "\n",
      "CPU times: user 2.34 s, sys: 8.12 ms, total: 2.35 s\n",
      "Wall time: 9.63 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "model = GradientBoostingRegressor()\n",
    "\n",
    "tuning_parameters = {\n",
    "    'learning_rate': [0.01, 0.05, 0.1],\n",
    "    'n_estimators' : [250, 500, 750, 1000, 1500],\n",
    "    'max_depth' : [2 ,3, 4],\n",
    "    'subsample' : [0.6, 0.8, 1.0]\n",
    "}\n",
    "\n",
    "# Using GridSearchCV would be too slow. Increase the number of iterations to explore more hyperparameter combinations.\n",
    "gb = RandomizedSearchCV(model, tuning_parameters, n_iter = 1, cv = 10, return_train_score=False, n_jobs=4)\n",
    "gb.fit(X_train, y_train)\n",
    "\n",
    "print('Best parameters found by randomised search:', gb.best_params_, '\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "GradientBoostingRegressor(alpha=0.9, ccp_alpha=0.0, criterion='friedman_mse',\n",
       "                          init=None, learning_rate=0.01, loss='ls', max_depth=2,\n",
       "                          max_features=None, max_leaf_nodes=None,\n",
       "                          min_impurity_decrease=0.0, min_impurity_split=None,\n",
       "                          min_samples_leaf=1, min_samples_split=2,\n",
       "                          min_weight_fraction_leaf=0.0, n_estimators=500,\n",
       "                          n_iter_no_change=None, presort='deprecated',\n",
       "                          random_state=None, subsample=1.0, tol=0.0001,\n",
       "                          validation_fraction=0.1, verbose=0, warm_start=False)"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gb.best_estimator_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our <TT>statlearning</TT> module contains a function for ploting the variable importances. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAm4AAAF1CAYAAABRUWbWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdebxe09n/8c8XJyQyKYJEOBJiJiQx10y11aIoqiktbfWn8VSLetrqozxF1aMloWqqak1FDE01iSGGmJKTyBypkpiqhgYVIjmS6/fHXnftnJw5x7nPvc/3/Xrdr7Pvtdde69p3/nBZa6+9FBGYmZmZWce3WrkDMDMzM7PmceJmZmZmViGcuJmZmZlVCCduZmZmZhXCiZuZmZlZhXDiZmZmZlYhnLiZWacm6VpJv25m3c0lhaQNG6kzUdLZqxDPeEnfb+31ZlZsa5Q7ADOzpki6B3g3Ir5Wz7kJwOyI+G5r2o6Ik1c1vrYUEQeXO4b6SPojsCgiTil3LGadmUfczKwS/BY4SlLvfKGkLYB9gKtb2qCk1SWpjeIrLP9OZh2LEzczqwRjgTeB4XXKvwU8FREzACT9QtJ8SYsk/V3SiFLF3DTn1yU9C3wArCvpj5KuytVrsI2cz0t6TtI7ku6StH5DgUuqljRa0j8l/UPSbySt3Uj9/0y15mL+mqS5kt6XNEZSb0kXS3pT0muSvp27/mRJz0r6Uerz9VR3jVydwZImSHpb0vOp7uqN/E7nAMcAJ6XfZVGqu5OkRyX9S9JCSX+RtFmunz9K+p2k69Jv9YqkFUY4Je0n6fEUy5uSrs2d20HS/ZLekvSSpJ9LqmrotzPrDJy4mVmHFxHLgWuBb5bKJHUBTmDF0bZZwB5AD+AU4JeSDqjT3HHA3kBP4O16umtOG18F9gQ2BVYHbqgvbkndgAnAdGAzYLv091eN3W89jgB2T/1tDkwC5gIbkv0mIyX1y9UfCGyQ+toL+BJweoppHeD+9NkA+AJZAnxanT7zv9PPgduA6yKie0R0T3WCLKnbCBgALAFurNPOMcBo4FPA94ErJW2cYtkJ+CtwVbqXTYE/pnMbAo+kfjci+zf5LHBmc380syJy4mZmleI6YGtJu6bvRwBVwJ9KFSLiDxHxWmQeIBupq5t0nRsRb0TEkohYVreTFrbxLnAW8DlJfeqJ+YtAbUT8LCIWR8RC4H+A4S2cfjwvIt6JiLeA+4DFEfG7iFgWEWOARcDgXP1a4OzU53PAJcDX07kvAO8DF0bE0oiYA/wSqPusX6O/E0BETIuIR1I77wDnAXtKWitX7f6I+EtELI+IP6W+d0znvgPclX7zJRHxQUQ8nM6dCEyOiGsjojYiXgF+Aaz0nKNZZ+LFCWZWESLiH5L+QjY69HT6+4eI+KBUR9L3gJOA/qmoG9kUa96CxvppRRul442BN+rU2wzYTNI7+S7Spw/wemOx5LyWO/6gzvdSWY/c99cjYnGdGDdOx/2B+RERufPP8/H95q9pVHrG8GJgl1z/AtYFXq0ndsgSt1LdauDJBprfDNinzm+3GrC8qbjMiswjbmZWSX4LHJOm2PYjN00qaR+yKb1vAutGRG+y0am6I1sN/oe/BW1U13P8Sj1NvgjMiYjeuU+viFgrIpqbtLXGBnVGvapz8b0MVNcZ8RuQyvPq/k71/W5Xk003bx8RPcmmVmHl36shC4AtGjj3IjC2zm/XM/2bmHVaTtzMrJKMA94C7gSejIhZuXM9gWVko2Mh6QtAS1+t0dw2fippfUm9gIvIEoy6o20A9wLdJf1QUndlNpZ0eAvjaqkq4EJJa0naHPgB8Pt07s9kI15nSeoiaWuy58aua6LNfwID6yR8Pcmmad9NCzR+1sI4rwK+JOkrKZZukvZN524Adpd0QrqP1SQNlPSZFvZhVihO3MysYqRFCteQTaPVfQXIfcAtwBSyxOtw4J4WdtHcNm4GngBeIhtdOrGBeBeRjQzuCMwD3iVbFLB9C+NqqefJ4l+Q4rwXuDTF9DZZMvo5sqnd+8iStsuaaPNqoDewMDd9+T1gf+DfwMNkSWGzRcRUsmfuTkvxvgh8JZ37R2r7qFS+kCxhr25JH2ZFoxUfczAzs0qWXrdxRkRsVe5YzKztecTNzMzMrEI4cTMzMzOrEJ4qNTMzM6sQHnEzMzMzqxBO3MzMzMwqhHdOKIBDDjkkxo4dW+4wzMzMrG00+BJrj7gVwFtvvVXuEMzMzKwdOHEzMzMzqxBO3MzMzMwqhF8HUgA9e68T/QZsW+4wzMzMOoX+ffswfszoT7KLBp9x8+KEAqhdWkvtsBHlDsPMzKxTeHnyyLL17anSdibpCEkhyfsImpmZWYs4cWt/xwETgWPLHYiZmZlVFidu7UhSd2BP4CRS4iZpNUlXSpotaYyk+yQdlc4NkfSIpCmSxknaqIzhm5mZWZk5cWtfhwNjI+JvwEJJOwNfAqqB7YGTgd0BJFUBI4GjImIIcD3w83IEbWZmZh2DFye0r+OAX6fjW9P3KuD2iFgO/FPShHR+S2A74H5JAKsDr7VvuGZmZtaROHFrJ5LWBfYHtpMUZIlYAHc1dAkwOyJ2b6cQzczMrIPzVGn7OQq4MSI2jYjqiOgPzAfeAo5Mz7ptAOyb6s8D1pf0n6lTSX5Zm5mZWSfmxK39HMfKo2t3An2BV4BZwG+Bp4F3I2IpWbL3C0nTgWnAHu0XrpmZmXU0niptJxGxbz1ll0O22jQiFqXp1EnAzHR+GrB3e8ZpZmZmHZcTt45hjKTeQBfg/Ij4Z0surupSRVUZ3+JsZmbWmfTv26dsfXuv0gIYOnRo1NTUlDsMMzMzaxsN7lXqZ9zMzMzMKoQTNzMzM7MK4cTNzMzMrEL4GbcC6Nl7neg3wK94M6sk/fv2YfyY0eUOw8w6pgafcfOq0gKoXVpL7bAR5Q7DzFrgZa8EN7NWqOipUknrSpqWPv+U9Grue5d66n9K0inNaHcNSe+k480lLU5tTpf0uKQt2iD2/SXtlvu+taRHUj9zJf0mlR8o6d3cfY1b1b7NzMysMlX0iFtE/AsYDCDpXGBRRFzSyCWfAk4BrmphV/MiotTPqcDZwEktDnhF+5Ntd/VU+j4KuDgi/qJsV/ntcnUnRMThq9ifmZmZVbiKHnFrjKSzJM1Kn9I84kXAlmnk6iJJPSU9JGmqpBmSDm1G0z2Bt1Mf20uanNqbIWlAGqGbJel6SbMl3SjpM5KekPQ3SUMlDQROBs5M1+4BbES29RWRmdn2v4qZmZlVsooecWuIpF2A44FdgNWBSZIeIRsp2zw3elYFHBYR70nqAzwOjKmnyS0lTSNL2tYEdk3l/w+4JCJuk7Qm2cOEGwNbAl8GngWmAksiYg9JRwJnR8RRkq4F3oqIX6dYLgUelfQ4MB74XUS8m/rZL/UPcGtEXNQmP5SZmZlVlKKOuH0auDMiPoiI94C7gb3qqSeyTdxnkCVL/SWtV0+9eRExOCIGAGfx8VTrE8BPJJ0F9I+ID1P53yNiTkQsB+YAD6TymUB1fQFHxLXANsAdwAHAk7nn9Cak/gc7aTMzM+u8ipq4NbiMto6vAb2AndMo3FvAWk1ccy9p4/eI+ANwBLAEuF9SaUP4Jbn6y3Pfl9PIKGdEvBoR10fEF8j+bbZu5n2YmZlZJ1DUxO1R4AhJXSV1Bw4DHgPeA3rk6vUC3oiIjyQdBPRrRtt7Ac8DSBoQEX+PiMuAvwA7tCDGFWKRdIikNdJxX2Ad4B8taM/MzMwKrpDPuEXEJEm3AJNT0W9KD/tLqpE0kyzRuhT4s6QasmfRnmugydIzbiIbPftWKv+KpOOAWrIk6ydAfVOt9bkHuF3Sl4BTgc8Cl0n6EAjgexHxZrbA1MzMzMw7JxRC127do9/w68odhpm1QNXkkcydOrHcYZhZx+SdE4qsqksVVX4Lu1lF6d+3T7lDMLMK5MStAAZtPpCaGv+fu5mZWdEVdXGCmZmZWeE4cTMzMzOrEE7czMzMzCqEV5UWQM/e60S/AduWOwwza6b+ffswfszocodhZh2XV5UWWe3SWmqHjSh3GGbWTC97FbiZtVKhpkolbSDpZkkvSJoi6UlJR5Q5pnskPVnOGMzMzKwYCpO4Kdti4G7g0YgYEBFDgGOBjZt5/eqfQEy9gZ2B3pI2a6CORz3NzMysWQqTuAH7A0sj4qpSQUS8GBEjJVVLekzS1PTZA0DSvpImSLoZKG2JdXcarZstqbS1FZJOkvQ3SQ9LukbSqFS+vqQ7JU1Onz1zMR0J/Bm4lSyJLLV1g6RLJU0AfiFpbUnXp+ufkXRYqldv3GZmZtY5FWm0Z1uy/Ubr8wZwUER8KGkL4BZgaDq3C7BdRMxP378REQsldQUmS7oTWBM4h2z07D3gIWB6qn8Z8KuImChpE2AcsHU6dxzwM+B14A7gwlxMg4ADI2KZpAuAhyLiG2mUbpKkB5qI28zMzDqZIiVuK5B0BbAXsBQ4EBglaTCwjCxpKpmUS9oATss9F9cf2ALYEHgkIhamtm/PtXEgsE1uM/ieknoA3YDNgYkREZI+krRdRMxK9W6PiGXp+GDgi5LOSN/XAjYh27i+objNzMyskylS4jabbGoSgIg4VdJ6QA1wOtmo145k08Mf5q57v3QgaV+yRGz3iPhA0sNkSVSDy3JTe7tHxOJ8oaSvA+sA81NS15NsuvQndftN7R8ZEfPqtHFuI3GbmZlZJ1OkZ9weAtaS9J1cWbf0txfwWkQsB4YDDS1E6AW8nZK2rYDdUvkkYB9J66TFBEfmrhkPfLf0JY2OQTZNekhEVEdENVBaLFGfccCItMACSTu1MG4zMzPrBAqTuEX2JuHDyRKs+ZImAb8HfghcCZwg6Smy6cb3G2hmLLCGpBnA+cBTqe1XgQuAp4EHgDnAu+ma04ChkmZImgOcIqmabKrzqVx884F/S9q1nn7PB6qAGZJmpe+0IG4zMzPrBLxzQjNJ6h4Ri9KI213A9RFxV7njAujarXv0G35ducMws2aqmjySuVMnljsMM+u4vHNCGzhX0oFkz7yNJ3tnXIdQ1aWKKr+J3axi9O/bp9whmFmF8ohbAQwdOjRqamrKHYaZmZm1jQZH3ArzjJuZmZlZ0TlxMzMzM6sQTtzMzMzMKoSfcSuAnr3XiX4Dti13GGaF079vH8aPGV3uMMys8/Gq0iKrXVpL7bAR5Q7DrHBe9mptM+tgCjdVKmkDSTdLekHSFElP5vYeLUc8n5VUI2mupGclXVKuWMzMzKyyFSpxS1tG3Q08GhEDIqK0zdTGzby+TbeUkrQdMAr4akRsDWwHvNCC6z0iamZmZv9RqMQN2B9YGhFXlQoi4sWIGCmpWtJjkqamzx6QbSwvaYKkm4GZqezuNFo3W9K3Sm1JOknS3yQ9LOkaSaNS+fqS7pQ0OX32TJecBfw8Ip5NsXwUEVema74g6WlJz0h6QNIGqfxcSVdLGg/cKGlbSZMkTUvbam3xif+KZmZm1iEVbURnW2BqA+feAA6KiA9T8nMLMDSd2wXYLu0nCvCNiFgoqSswWdKdwJrAOcDOwHtkm9pPT/UvA34VERMlbUK2aXxphO3/GohnIrBbRISkk8mSvB+kc0OAvSJisaSRwGURcZOkLnijeTMzs06raInbCiRdAewFLAUOBEZJGgwsI9u0vWRSLmkDOC33XFx/YAtgQ+CRiFiY2r4918aBwDbZTC0APSX1aCK8jYHbJG0EdAHy/d8bEYvT8ZPAjyVtDIyOiOeacetmZmZWQEWbKp1NNiIGQEScChwArA+cDrwO7Eg20tYld937pQNJ+5IlYrtHxI7AM2T7kza4NJfsd9w9IganT7+IeC/FM6SBa0YCoyJie+DbqY+V4omIm4EvAouBcZL2byQOMzMzK7CiJW4PAWtJ+k6urFv62wt4LSKWA8NpeMqxF/B2RHwgaStgt1Q+CdhH0jpp0cCRuWvGA98tfUmjegC/BH4kaVAqX03S93P9vJqOT2johiQNAF6IiMuBe4EdGqprZmZmxVaoxC2ytwkfTpZgzZc0Cfg98EPgSuAESU+RTXG+30AzY4E1JM0AzgeeSm2/ClwAPA08AMwB3k3XnAYMTYsH5gCnpGtmAN8DbpE0F5gFbJSuORe4XdJjwFuN3NYxwCxJ04CtgBub/4uYmZlZkXjnhBaQ1D0iFqURt7uA6yPirnLH1bVb9+g3/Lpyh2FWOFWTRzJ36sRyh2FmnY93Tmgj50o6kOx5tPFk74wru6ouVVT5De9mba5/3z7lDsHMbAUecSuAoUOHRk1NTbnDMDMzs7bR4IhboZ5xMzMzMysyJ25mZmZmFcKJm5mZmVmF8DNuBdCz9zrRb8C25Q7DrEPq37cP48eMLncYZmYt4VWlRVa7tJbaYSPKHYZZh/SyV1ybWYF4qrQRkq6X9IakWU3U21fSHrnv50p6VdK09LkolT8saWgDbRwq6RlJ0yXNkfTtxtoyMzOzzscjbo27ARhF07sV7AssAp7Ilf0qIi5pTieS1gSuBnaJiFfS9+rWtGVmZmbF5RG3RkTEo8DCfJmk09KI2AxJt0qqJtvi6vQ0Ivbp5rQtaZGk8yQ9DexKlkT/K/W7JCLmteW9mJmZWeVz4tZyZwM7RcQOwCkRsQC4imxUbHBEPJbqnZ6b3vxMPe2sDcyKiF1Tgngv8KKkWyQdLyn/b9NUW2ZmZtYJOHFruRnATZK+CnzUSL1SIjc4IsbVc34ZcGfpS0ScDBwATALOAK5vQVtmZmbWCThxa7nPA1cAQ4ApacP51vgwIpblCyJiZkT8CjgIOHLVwjQzM7OiceLWAmn6sn9ETADOAnoD3YH3gB6r0G53SfvmigYDL65CqGZmZlZAXlXaCEm3kK0YXU/SK8D5wHBJvchejveriHhH0p+BOyQdBrTmhWoCzpL0W2Ax8D5wYhvcgpmZmRWId04ogK7duke/4deVOwyzDqlq8kjmTp1Y7jDMzFqiwZ0TPFVqZmZmViE8VVoAVV2qqPK2Pmb16t+3T7lDMDNrM07cCmDQ5gOpqfFUkJmZWdF5qtTMzMysQjhxMzMzM6sQXlVaAD17rxP9Bmxb7jDMyqJ/3z6MHzO63GGYmbWlBleV+hm3AqhdWkvtsNa8Ps6s8r3shTlm1ol4qtTMzMysQjhxa4CkZZKm5T7VTdRfIGm9dLwo/a2WtDhdP13SE5K2bKKdaklfyX0/UdKoVb8jMzMzq3RO3Bq2OCIG5z4LWtnO8+n6HYHfAz9qon418JUm6piZmVkn5MStBeqOfkkaU2dz+Kb0BN5O11ZLekzS1PTZI9W5CPh0GqU7PZX1lTRW0nOSLm6LezEzM7PK48UJDesqaVo6nh8RR7SynYGpnR5AN2DXVP4GcFBEfChpC+AWYChwNnBGRBwKWbIIDAZ2ApYA8ySNjIiXWxmPmZmZVSgnbg1bHBGD26Cd50vtSDoGuBo4BKgCRkkaDCwDBjXSxoMR8W5qYw6wKeDEzczMrJNx4tYyH7Hi9PJaLbz+XuB36fh04HVgx9Tmh41ctyR3vAz/u5mZmXVKfsatZRYAgyWtJqk/sEsLr98LeD4d9wJei4jlwHBg9VT+Htm0qpmZmdkKPHLTMo8D84GZwCxgajOuKT3jJmApcHIqvxK4U9LRwATg/VQ+A/hI0nTgBtJiBjMzMzNveVUAXbt1j37Dryt3GGZlUTV5JHOnTix3GGZmbclbXhVZVZcqqrztj3VS/fv2KXcIZmbtxolbAQzafCA1NR5xMDMzKzovTjAzMzOrEE7czMzMzCqEEzczMzOzCuFVpQXQs/c60W/AtuUOw2wl/fv2YfyY0eUOw8ys0nhVaZHVLq2ldtiIcodhtpKXvdrZzKxNeao0R1J/SRMkzZU0W9J/tfD6hyUNTccLJM2UNC199pBULWlWA9euJulySbPSdZMlbdZQW6t+t2ZmZlZpPOK2oo+AH0TEVEk9gCmS7o+IOa1sb7+IeKv0RVJ1fZUkrQEcDfQFdoiI5ZI25uPdFFZqy8zMzDofJ245EfEa8Fo6fk/SXKCfpCuBp4H9gN7ASRHxmKSuZJvGbwPMBbo2ty9JJwKfJ9uofm1gDB/vXUpEvNJW92VmZmbF4MStAWl0bCeyhA1gjYjYRdLngP8BDgS+A3wQETtI2oGV9y6dIGkZsCQidq2nm93JRtgWphG2iZI+DTwI/DEinmlBW2ZmZlZwTtzqIak7cCfwvYj4tySA0tK4KUB1Ot4buBwgImZImlGnqaamN++PiIXp+lckbQnsnz4PSjo6Ih5sZltmZmZWcE7c6pBURZa03RQR+fcYLEl/l7Hi77Yq71PJP8NGRCwB/gr8VdLrwOFko29mZmZmXlWap2xo7TpgbkRc2oxLHgWOT9duB+ywCn3vLKlvOl4ttfVia9szMzOz4vGI24r2BIYDMyVNS2U/aqT+b4DfpSnSacCkVei7D3CNpDXT90nAqFVoz8zMzArGOycUQNdu3aPf8OvKHYbZSqomj2Tu1InlDsPMrNJ454Qiq+pSRZXfUG8dUP++fcodgplZoThxK4BBmw+kpsajGmZmZkXnxQlmZmZmFcKJm5mZmVmFcOJmZmZmViG8qrQAevZeJ/oN2LbcYVgF6N+3D+PHjG66opmZlZNXlRZZ7dJaaoeNKHcYVgFe9upjM7OK1ummSiUtkzQt96mWNFTS5W3YxwJJ67VVe2ZmZmbQOUfcFkfE4DplC4CauhUlrRERH7VLVGZmZmZN6HQjbvWRtK+kMen4XElXSxoP3ChpdUm/lDRZ0gxJ385d86ikuyTNkXRV2mO0btt3S5oiabakb+XKD5E0VdJ0SQ+msrUlXZ/6ekbSYal8W0mT0gjhDElbtMsPY2ZmZh1KZxxx65rbh3R+RBxRT50hwF4RsTglW+9GxLC0j+jjKakD2AXYhmwz+LHAl4A76rT1jYhYKKkrMFnSnWQJ8zXA3hExX9KnUt0fAw9FxDck9QYmSXoAOAW4LCJuktQFWL0tfggzMzOrLJ0xcatvqrSueyNicTo+GNhB0lHpey9gC2ApMCkiXgCQdAuwFysnbqdJKiWH/dO16wOPRsR8gIhYmOvri5LOSN/XAjYBngR+LGljYHREPNeiOzYzM7NC6IyJW3O8nzsWMCIixuUrSNoXqPsulainzoHA7hHxgaSHyZIx1XNtqa8jI2JenfK5kp4GPg+Mk3RyRDzUojsyMzOziudn3Jo2DviOpCoASYMkrZ3O7SJps/Rs2zFA3Q1DewFvp6RtK2C3VP4ksI+kzVKbpanSccAISUrlO6W/A4AXIuJy4F5gh0/iRs3MzKxjc+LWtGuBOcBUSbOA3/LxSOWTwEXALGA+cFeda8cCa0iaAZwPPAUQEW8C3wJGS5oO3Jbqnw9UATNSX+en8mOAWenZvK2AG9v6Js3MzKzj884JrZSmQc+IiEPLHUvXbt2j3/Dryh2GVYCqySOZO7XuwLCZmXUw3jmhyKq6VFHlN+JbM/Tv26fcIZiZ2SrwiFsBDB06NGpqVnp/sJmZmVWmBkfc/IybmZmZWYVw4mZmZmZWIZy4mZmZmVUIP+NWAD17rxP9Bmxb7jCsA+rftw/jx4wudxhmZtYyXlVaZLVLa6kdNqLcYVgH9LJXG5uZFYqnSs3MzMwqREUlbpKWSZqW+5zdRP0ftbKfLpJ+Lel5SX+XNEbSJq2LGiSdm9s4vu65K9K9zJG0OHdvR6Xza0h6S9KFre3fzMzMiqHSpkoXR8TgFtT/EXBBSzqQtHq6pgcwKCKWSfo6cI+kIRGxvCXtNSUiTk39VgNj6rm/g4F5wJcl/Sj8UKKZmVmnVVEjbvWR1EvSPElbpu+3SPqmpIuArmn06qZ07quSJqWy36YkDUmLJJ0n6WlgT+DrwOkRsQwgIn4HLAIOlFSd9hEt9X+GpHPT8TclTZY0XdKdkrq1wS0eB1wGvMTHm9SbmZlZJ1RpiVspESt9jomId4HvAjdIOhZYJyKuiYizSSN0EXG8pK3JNmvfM41qLQOOT+2uDcyKiF2Bd4CXIuLfdfquAbZpIr7RETEsInYE5gInrcrNSuoKHACMAW4hS+LMzMyskyrEVGlE3C/paOAKYMcGrj0AGAJMlgTQFXgjnVsG3JmOBdQ3Hdng0tyc7ST9L9Ab6A6Ma8Y1jTkUmBARH0i6EzhH0n9GAs3MzKxzqbTErV6SVgO2BhYDnwJeqa8a8PuI+O96zn2YS4b+DmwqqUdEvJerszNwB/ARK45UrpU7vgE4PCKmSzoR2Lfld7OC44A9JS1I39cF9gMeWMV2zczMrAJV2lRpQ04nm5o8DrheUlUqr80dPwgcJakPgKRPSdq0bkMR8T7we+DS3DNwXwM+BB4HXgf6SFpX0ppko2IlPYDXUp/Hswok9QT2AjaJiOqIqAZOxdOlZmZmnValjbh1lTQt930scD1wMrBLRLwn6VHgJ8D/AFcDMyRNTc+5/QQYn0boaskSoRfr6ee/gV8C89JzZm8Cu6cVnbWSzgOeBuYDz+auOyeVvwjMJEvkWutLwEMRsSRXdg9wsaQ165SbmZlZJ+Atr5ogaUOyBPHKiLi63PHUx1teWUO85ZWZWUVq8Ll6J24FMHTo0KipqSl3GGZmZtY2vFdpRyDpCrL3xOVdlt4TZ2ZmZtYoJ27tqLRLgpmZmVlrFGVVqZmZmVnh+Rm3AvDihI7PiwTMzKwF/IxbkdUuraV22Ihyh2GNeHnyyHKHYGZmBeCpUjMzM7MK0WTiJmlZnY3dqyUNlXR5WwUhaYGk9dqqvdTmXpImSXpW0jxJrVoYkO43JJ2fK1tPUq2kUen7KWl3hZa021fSHa2JyczMzDqn5kyV1rex+wJgpReHSVojIj5qi8BWRXpp7s1k+4ZOTUnhOEn/iIi7WtHkC2RbW52Tvh8NzC6djIirWtpgRPwDOKoVsZiZmVkn1aqpUkn7ShqTjs+VdLWk8cCNklaX9EtJkyXNkPTt3DWPSrpL0hxJV6Wtp+q2fbekKZJmS/pWrvwQSVMlTZf0YCpbW9L1qa9nJB2Wqp8K3BARUwEi4i3gLODMdN0Nko7Ktb2oiVteDMyVNDR9Pwb4U+76cyWdkY5PS/c3Q9KtqWyf3IjlM5J6pJG8WaEVa+4AACAASURBVOn8iZJGSxor6TlJF+faPknS3yQ9LOma0iifmZmZdT7NGXHL7w86PyKOqKfOEGCviFickq13I2JY2oT98ZTUAewCbEO2l+dYsv04604XfiMiFqY9QidLupMswbwG2Dsi5kv6VKr7Y7L9PL8hqTcwSdIDwLZkG8Xn1aS+W+tW4FhJ/wSWAf8A+tZT72xgs4hYkmICOAM4NSIel9SdbMP6ugYDOwFLyPZIHZn6OQfYGXgPeAiYvgr3YGZmZhWstVOldd0bEYvT8cHADrkRrV7AFsBSYFJEvAAg6RZgL1ZO3E6TVEoO+6dr1wcejYj5ABGxMNfXF0ujXcBawCZky2jb+j0nY4HzgdeB2xqpNwO4SdLdwN2p7HHgUkk3AaMj4hVppZW+D0bEuwCS5gCbAusBj5TuV9LtwKA2uh8zMzOrMG21qvT93LGAERExOH02i4jSiFvdZGqF75L2BQ4Edo+IHYFnyJKxhhIxAUfm+tokIuaSPX82tE7dIXz8XN5HpHtXlkF1aeoGI2IpMAX4AXBnI1U/D1yR+puSnvu7CDgZ6Ao8JWmreq5bkjteRpZUN/geFzMzM+t8PonXgYwDviOpCkDSIElrp3O7SNosPdt2DDCxzrW9gLcj4oOU3OyWyp8E9pG0WWqzNFU6DhiRki8k7ZTKrwBOlDQ4la8L/JxsxAyyxRVD0vFhQFUz7+3/gB9GxL/qO5nuq39ETCB7pq430F3SwIiYGRG/IEse60vc6jOJ7L7XkbQGcGQzrzMzM7MC+iRewHstUA1MTQnVm8Dh6dyTwEXA9sCjQN0VnmOBUyTNAOYBTwFExJvp2bnRKTl6AziILBH7NTAj9bUAODQiXpP0VeBqSb1SPCdGxCOpn2uAeyRNAh5kxRHDBkXEbHKrSeuxOvDH1KeAX0XEO5LOl7Qf2UjaHOCvwEbN6O9VSRcAT5M9UzcHeLc5sZqZmVnxtNuWV2ka9IyIOLRdOlyx71OBU8gWN7zd3v2vCkndI2JRGnG7C7i+7itNvOVVx+ctr8zMrAU695ZXEXEF2fRpJTpX0oFkz/qN5+MFD/8xaPOB1NTUnXU2MzOzomm3xC0iHgYebq/+WkrS9sAf6hQviYhdyxFPSUSc0XQtMzMz6ww6xYhbc0TETLJ3qZmZmZl1SN5k3szMzKxCtNviBPvkdPbFCX7w38zMCqZzL04outqltdQOG1HuMMrm5ckjyx2CmZlZu/BUqZmZmVmFKHviJmlR7vhzkp6TtImkUyR9LZWfKKm+Dd3z7ZwoaVQbxnW4pBmSnpU0K7f3amvaqpY0q4Fzn5E0LX0WSZqXjm/M1blM0qvp5cNmZmbWSXWYqVJJBwAjgYMj4iXgqtzpE4FZZLsHtEcsOwKXAAdFxPy01dYDkuZHxJS27CsixpFt3YWkh8leUlzaU7W0jdYRwMvA3nTgV6qYmZnZJ6tDjOBI+jTZNlSfj4jnU9m5ks5II11DgZvSSFRXScMkPSFpuqRJknqkpvpKGptG7S7OtX+wpCclTZV0u6TuqXyBpJ+l8pm5zd/PAC6IiPkA6e8FZBvMI+lhSUPT8XqSFqTjakmPpfamStqjDX6e/ciS1t8Ax7VBe2ZmZlahOkLitiZwD3B4RDxb92RE3EG2MfvxETGYbL/P24D/iogdgQOBxan6YLLN67cHjpHUX9J6wE+AAyNi59TW93NdvJXKf0OWsAFsC9QdWasBtmniXt4gG6XbOcVxeVM33wzHAbeQbXd1qKSqNmjTzMzMKlBHSNxqgSeAk5pZf0vgtYiYDBAR/46Ij9K5ByPi3Yj4kGxD9k2B3cgSrsclTQNOSOUlpfdITCHbjB6yZbh135PS4NLcnCrgGkkzgdtpOtFrlKQuwOeAuyPi32SbzR+8Km2amZlZ5eoIz7gtB75M9gzZjyLigibq15dUlSzJHS8juz8B90dEQ9OMS+rUB5hNNj07I1evNFoH8BEfJ71r5eqcDrwO7JjOf9jYjTTDIUAvYKYkgG7AB8BfVrFdMzMzq0AdYcSNiPgAOBQ4XlJ9I2/vAaXn2J4le5ZtGICkHpIaS0CfAvaUtHmq303SoCZCugT4b0nV6Zpq4HvAL9P5BcCQdJxfbdqLbDRwOTAcWL2JfppyHHByRFRHRDWwGXCwpG6r2K6ZmZlVoA6RuAFExEKyEaafSDqszukbgKvSVOfqZM+PjZQ0HbifFUe96rb7Jtmq1FskzSBL5LZqqH66ZhrwQ+DPkv4G/A34TkTMS1UuAb4j6QlgvdylVwInSHoKGAS839R9NyQlZ58hN7oWEe8DE4EvtLZdMzMzq1ze8qoZJF0E7Ap8JiKWljueurzllbe8MjOzQmnwuXonbgUwdOjQqKmpabqimZmZVQLvVdoRSPoM8Is6xfMj4ohyxGNmZmaVxYlbO8rvkmBmZmbWUh1mcYKZmZmZNc7PuBVAZ1yc4AUJZmZWYH7Grchql9ZSO2xEucNoVy9PHlnuEMzMzNqdp0qbIGlRC+oeLmmbOmVrSHpL0oVtH52ZmZl1Jk7c2tbhrLw/6cHAPODLSvtW1SVpVXdYMDMzs07AiVsrSNpU0oOSZqS/m0jaA/gi8EtJ0yQNTNWPAy4DXiLb8L7UxgJJP5U0ETha0kBJYyVNkfSYpK1SvS9IelrSM5IekLRBO9+umZmZdRBO3FpnFHBjROwA3ARcHhFPAPcCZ0bE4Ih4XlJX4ABgDHALWRKX92FE7BURtwJXAyMiYghwBtn2WZBtcbVbROwE3Aqc9UnfnJmZmXVMXpzQOrsDX0rHfwAubqDeocCEiPhA0p3AOZJOj4hl6fxtAJK6A3sAt+dmU9dMfzcGbpO0EdAFmN+md2JmZmYVw4lb22jonSrHAXtKWpC+rwvsBzyQvpc2oV8NeCciBtfTxkjg0oi4V9K+wLltEbCZmZlVHk+Vts4TwLHp+Hiy6UyA94AeAJJ6AnsBm0REdURUA6ey8nQpEfFvYL6ko9O1krRjOt0LeDUdn9D2t2JmZmaVwolb07pJeiX3+T5wGvB1STOA4cB/pbq3AmdKegY4GngoIpbk2roH+KKkNVnZ8cBJkqYDs4HDUvm5ZFOojwFvtfXNmZmZWeXwzgkF0LVb9+g3/Lpyh9GuqiaPZO7UiU1XNDMzqzwN7pzgETczMzOzCuHFCQVQ1aWKqk62BVT/vn3KHYKZmVm7c+JWAIM2H0hNjacNzczMis5TpWZmZmYVwombmZmZWYXwqtIC6Nl7neg3YNtyh9Fu+vftw/gxo8sdhpmZ2SelwVWlfsatAGqX1lI7bES5w2g3L3eyhRhmZmYlnio1MzMzqxBO3MzMzMwqRGESN0kbSrpV0vOS5ki6T9KgVrRzoqS+rbjuXEln5L6vIektSRfWqXetpG2a2eYFkn6R+76ppBck9W5pfGZmZlb5CpG4SRJwF/BwRAyMiG2AHwEbtKK5E4F6EzdJq7egnYOBecCXU3wARMTJETGnmW2fDxwmaev0/TLgnIh4pwVxmJmZWUEUInED9gNqI+KqUkFETIuIxySdKWmypBmSfgYgqVrSXEnXSJotabykrpKOAoYCN0malsoWSPqppInA0ZK+mdqbLulOSd0aiOk4skTrJWC3UqGkhyUNTceLJJ0n6Wlg97oNRMRi4PvAlZI+C/SIiJva4gczMzOzylOUxG07YErdQkkHA1sAuwCDgSGS9k6ntwCuiIhtgXeAIyPiDqAGOD4iBqfECeDDiNgrIm4FRkfEsIjYEZgLnFRPv12BA4AxwC1kSVx91gZmRcSuEVHv1gcRcR+wELgR+H9N/RBmZmZWXEVJ3BpycPo8A0wFtiJL2ADmR8S0dDwFqG6kndtyx9tJekzSTOB4oL4XqB0KTIiID4A7gSMamApdls435QpgckTMa0ZdMzMzK6iivMdtNnBUPeUCLoyI365QKFUDS3JFy4CujbT/fu74BuDwiJgu6URg33rqHwfsKWlB+r4u2XTuA3XqfRgRyxrpt2R5+piZmVknVpQRt4eANSV9s1QgaRjwb+Abkrqnsn6S+jTR1ntAj0bO9wBek1RFNuK2Akk9gb2ATSKiOiKqgVNpeLrUzMzMrFkKMeIWESHpCODXks4GPgQWAN8je37tybSwcxHwVbIRtobcAFwlaTH1LBgAzgGeBl4EZrJykvcl4KGIyI/o3QNcLGnNlt2ZmZmZ2ce8V2kBeK9SMzOzQvFepUU2aPOB1NTUuyjVzMzMCsSJWwch6S5gszrFP4yIceWIx8zMzDoeJ24dREQcUe4YzMzMrGMryqpSMzMzs8Lz4oQC6AyLE7wgwczMOhEvTiiy2qW11A4bUe4wPlEvTx5Z7hDMzMzKrnBTpZKWpQ3ip0uaKmmPNmhzsKTP5b6fKOnN1M80STem8vMkHdhEWxtIGpPimyPpvlReLWlxrs1pkrpI2krSk5KWSDpjVe/FzMzMKlcRR9wWR8RgAEmfAS4E9lnFNgcDQ4H7cmW3RcR385Ui4qfNaOs84P6IuCzFuEPu3POl2EskLQROAw5vTeBmZmZWHIUbcaujJ/A2gKSNJD2aRrJmSfp0Kl8k6ReSpkh6QNIukh6W9IKkL0rqQpZsHZOuPaahziTdIOmodLxA0s/SqN9MSVulahsBr5SuiYgZjd1ARLwREZOB2lX5IczMzKzyFTFx65oSrGeBa4HzU/lXgHFpRGtHYFoqXxt4OCKGkO1T+r/AQcARwHkRsRT4KdkI2+CIuC1dV0rkpkn6egOxvBUROwO/AUrTnFcA10maIOnHkvrm6g/MtXnFqv4QZmZmVixFnyrdHbhR0nbAZOD6tDn83RFRStyWAmPT8UxgSUTUSpoJVDfSz0pTpfUoLYOcQraHKRExTtIA4BDgs8AzKT6oZ6rUzMzMrKSII27/ERFPAusB60fEo8DewKvAHyR9LVWrjY/fibIcWJKuXc6qJ7aljeaX5duKiIURcXNEDCdLKPdexX7MzMysEyh04paeK1sd+JekTYE3IuIa4Dpg5xY09R7Qo41i2l9St3TcAxgIvNQWbZuZmVmxFXGqtKuk0jSogBMiYpmkfYEzJdUCi4CvNdRAPSYAZ6d2L1zF+IYAoyR9RJY4XxsRkyVV11dZ0oZADdlCi+WSvgdsExH/XsU4zMzMrMJ454QC6Nqte/Qbfl25w/hEVU0eydypE8sdhpmZWXvwzglFVtWliqqC7yzQv2+fcodgZmZWdk7cCmDQ5gOpqfFolJmZWdEVenGCmZmZWZE4cTMzMzOrEE7czMzMzCqEV5UWQM/e60S/AduWO4xm69+3D+PHjG66opmZWefkVaVFVru0ltphI8odRrO9XPAVsGZmZp8UT5WamZmZVYgOkbhJCkl/yH1fQ9Kbksak7xtIGiNpuqQ5ku5L5atJulzSLEkzJU2WtFkTfd0g6agGzu0i6VFJ8yQ9K+laSd0knShpVFvecz19j033N1vSVZJWT+WfknS/pOfS33U+yTjMzMys4+oQiRvwPrCdpK7p+0Fkm8GXnAfcHxE7RsQ2wNmp/BigL7BDRGwPHAG805oAJG0A3A78MCK2BLYGxtJGe5Q2w5cjYkdgO2B94OhUfjbwYERsATzIx/duZmZmnUxHSdwA/gp8Ph0fB9ySO7cR8ErpS0TMyJW/FhHLU/krEfE2gKRFpfqSjpJ0Q669AyU9Julvkg5NZacCv4+IJ1NbERF3RMTr+SAlfUHS05KekfRASviQtI+kaenzjKQekjZKI3jT0qjgpxu6+dzeo2sAXYDSqpHDgN+n498DhzfUhpmZmRVbR0rcbgWOlbQWsAPwdO7cFcB1kiZI+rGkvqn8T8AXUmL0f5J2amZf1cA+ZIniVanP7YApzbh2IrBbROyUYj4rlZ8BnBoRg4FPA4uBrwDjUtmOwLTGGpY0DngDeA+4IxVvEBGvAaS/3vvJzMysk+owiVsaRasmG227r865ccAA4BpgK+AZSetHxCvAlsB/A8uBByUd0Izu/hQRyyPiOeCF1GZzbQyMkzQTOBMovYfjceBSSacBvSPiI2Ay8HVJ5wLbR8R7jTUcEZ8hG0VcE9i/BTGZmZlZJ9BhErfkXuASVpwmBSAiFkbEzRExnCwh2juVL4mIv0bEmcAFfDyVmH9B3Vp1m6vn+2xgSDNiHAmMSs/UfbvUdkRcBJwMdAWekrRVRDya4nwV+IOkrzXVeER8SPY7HJaKXpe0EUD6+0YzYjQzM7MC6miJ2/XAeRExM18oaX9J3dJxD2Ag8JKknUvTppJWI5tifTFd9rqkrVP5EXX6OTqtSB1INpI3DxgFnCBp11y/X5W0YZ1re/HxwokTcnUHRsTMiPgFUANsJWlT4I2IuAa4Dti5vpuW1D2XnK0BfA54Np2+N9fPCcA99bVhZmZmxdehXsCbpj4vq+fUEGCUpI/Iks1rI2KypEOAayStmepNIkvAIFt9OQZ4GZgFdM+1Nw94BNgAOCWNcn0o6VjgEkl9yKZeHwXqvuL/XOB2Sa8CTwGl1498T9J+wDJgDtlii2OBMyXVAouAhkbc1gbuTfexOvAQcFU6dxHwJ0knAS/x8WpTMzMz62S85VUBeMsrMzOzQvGWV0U2aPOB1NRMLHcYZmZm9glz4tbOJD1Ntmo0b3jd5/rMzMzM6nLi1s4iYtema5mZmZmtrKOtKjUzMzOzBnhxQgG09+IELy4wMzP7RHlxQpHVLq2ldtiIduvv5ckj260vMzMz+5inSs3MzMwqRFkTN0kbSLpZ0guSpkh6UlLdXQ7aM57PSqqRNFfSs5IuaaN2b5B0VAPn7pI0TdLfJb2bjqdJ2iOdX19SraRvt0UsZmZmVrnKlrhJEnA38GhEDIiIIWQ7DWzczOtXb+N4tiPbdeGrEbE1sB3ZBvSfqIg4IiIGk+1z+lhEDE6fJ1KVo8l2aDjuk47FzMzMOrZyjrjtDyyNiNLWTkTEixExUlK1pMckTU2f0ujTvpImSLoZmJnK7k6jdbMlfavUlqSTJP1N0sOSrpE0KpWvL+lOSZPTZ890yVnAzyPi2RTLRxFxZbpmU0kPSpqR/m6Sym+QdLmkJ9Ko4VGpXJJGSZoj6S9An1X4nY4DfgBsLKnfKrRjZmZmFa6cidu2wNQGzr0BHBQROwPHAJfnzu0C/Dgitknfv5FG64YCp0laN208fw6wG3AQsFXu+suAX0XEMOBI4NpUvh0wpYF4RgE3RsQOwE114tkI2As4lGxfUcg2td8S2B74JrBHA+02SlJ/YMOImAT8iey3MDMzs06qw6wqlXQFWQK0FDiQbFP5wWSbtg/KVZ0UEfNz30/LPRfXH9gC2BB4JCIWprZvz7VxILBNNlMLQE9JPZoIb3fgS+n4D8DFuXN3R8RyYI6kDVLZ3sAtEbEM+Iekh5povyHHkiVsALcC1wGXtrItMzMzq3DlTNxmk414ARARp0paD6gBTgdeB3YkGxX8MHfd+6UDSfuSJWK7R8QHkh4G1qKR95+k9naPiMX5QkmzgSHA9GbEnn/53ZJ8Mw3Uaa3jgA0kHZ++95W0RUQ81wZtm5mZWYUp51TpQ8Bakr6TK+uW/vYCXksjWcOBhhYi9ALeTknbVmRTowCTgH0krSNpDXIJIjAe+G7pSxrVA/gl8CNJg1L5apK+n849QTb6BXA80NSO7o8Cx0paXdJGwH5N1F+JpC2BtSOiX0RUR0Q1cGEuDjMzM+tkypa4/f/27jxIr6pO4/j3mZBIIAYkEjQLhiWAoJCQgGQEBaRcGDVQ7AMCMpYyIiAziBGnFCcFhdsoAsMUCoKK7IuAU4BsAs4A6SyGJcSgIAExyL5DA8/8cU/Dpac73Ql5++337edT1ZX3nnvOuee9uZ386iz3uNqyYTeqAOs+SbcDZwNfBf4TOEjSrVRDnM/1Us1VwGqSFgKzqVZfYvsh4ATgNuBa4G7gqVLmCGB6WWhwN3BoKbMQ+DJwrqRFwJ1U89e6yny2XOczwJF9fL1LgSVUCyhOA37br5vyZvuVeuouJqtLIyIihqy23fJK0ijbz5Yet0uBM213D4TaQra8ioiIaCtDcsur4yTtQjXn7Rqqd8a1pU023oiOjr5GbyMiIqLVtW3gZvvoZrehO0mXAht0S/6q7aub0Z6IiIhoLW0buA1Gtpu2nVdERES0vmwyHxEREdEi2nZxwlAykIsTsjAhIiKi4Ybk4oQho/PlTjq3OXxArrV0zskDcp2IiIj4/zJUGhEREdEiWiJwkzRB0q8kLZH0R0knSRrR4Gs+W/6cJOnOWvr2km6XdI+kxZIOWxXX6SPPaEkPSTrlrVwrIiIiWtugD9xU7QZ/CdVm7pOpdlIYBRz/Futd4WFiSe8Cfgkcansz4IPAIbVN7htlNiu3+0JERES0kUEfuAE7Ay/a/imA7VepNqE/RNIcSa/Pypd0o6RpktaUdGY5P1/SzHL+YEkXSroCuEbSKEnXSZon6Y6ufMtxGHCW7XmlLY8CxwBfKfWfJWnPWnu6eu1W9DqvkzQNWI/qJcIRERExhLVC4LYFMLeeYPtp4AHgSmBvgLKZ+zjbc4GvA9fb3oZqg/fvSlqzFJ8BHGR7Z+BFYHfbW5d83y89fP1uC9ABbN7Hd1jR61C+098B36cEhhERETG0tULgJqCnd5YIuBHYqxzvDVxYPn8UmCVpQcmzOrB+Ofcb24/X6jihbB5/LTCeqndrRdvSn++wItfp8kXgv20vXYlrRkRERJtphdeB3AXsUU+QNBqYCMwBHpO0JbAP8IWuLMAethd3K/cB4Lla0v7AusA0252S7qcK8pbXlunA5bW0aVS9bgCvUILh0qPWtYBiRa/TZQawg6QvUs3rGyHpWduz+lE2IiIi2kwr9LhdB6wh6UAAScOohg/Psv08cB7VPLO1bN9RylwNHN41HClpai91rwU8UoKpnYD39NGWU4GDJU0p9Y6hWiQxu5y/nyqQA5gJDF/J6wBge3/b69ueBBwN/CxBW0RExNA16AM3V1s77A7sJWkJ8AeqOWPHliwXAfsCF9SKzaYKmhaWV3nMpmfnANMldVD1it3TR1seBg4ATpe0GPgL8CPbXSs+fwx8WNLtQL13b4WuExEREdGTbHn1FpR3uB0KfMj2E81qR7a8ioiIaCu9LmBM4NYGpk+f7o6Ojr4zRkRERCvIXqWDmaT3Az/vlvyS7Q80oz0RERExOCVwGwTKooopzW5HREREDG6DfnFCRERERFQyx60NNHpxQhYkREREDKjMcWtnnS930rnN4Q2rf+mckxtWd0RERPRfhkprujaFrx0fLOmUPsq8nkfSupJuKxvb7yDp/rKp/IL+bi4v6dja50nlPXQRERERCdxWsY8A99ieavvmkraT7SnAnsCP+lHHsX1niYiIiKEogVs/SfpUrTftWknrdTs/BfgOsGvpYRvZrYrRwBO1/JdJmivpLkmfL2knAiNL+XNK1mGSflzyXdNDvRERETFEJHB7s66gaYGkBcC/187dAmxneypv7I/6OtsLgG8A59ueYvuFcuqGMtz5W+DfakUOsT2NatP6IySNKfuQvlDK71/yTQZOtb0F8CSwx6r9yhEREdEqsjjhzV4ow5pANX+NKrACmACcL+ndwAjgvn7WuZPtRyVtBFwn6Ubbz1IFa7uXPBOpArTHeih/XwkKAeYCk1bkC0VERET7SI9b/50MnGL7/cAXgNVXpLDtPwLLgM0l7QjsAsywvRUwfzn1vVT7/CoJtiMiIoasBG79txbwUPl80IoWljQW2AD4c6nrCdvPS9oM2K6WtVPS8Lfa2IiIiGg/Cdz67zjgQkk3A4+uQLkbyny5G4BZtpcBVwGrSVoIzAZureU/HVhYW5wQERERAWTnhLYwco1RHv+ZMxpW//A5J7No3i0Nqz8iIiLeJDsntLPhI4YzvIG7G0wcN7ZhdUdERET/JXBrA5tsvBEdHekRi4iIaHeZ4xYRERHRIhK4RURERLSIBG4RERERLSKrStvA6LXf4fEbbrHK6504bizXXHnJKq83IiIiliurSttZ58uddG5z+Cqvd2kDV6pGRETEistQaURERESLGFSBm6QxkhaUn79Keqh2PKKH/OtIOrR2vLGkF0r+RZLOkrTKehUl/brsnFBP+4Wk3Vawnl0lzZF0T2nruZIm9KPcapKeXNF2R0RERHsYVIGb7cdsT7E9Bfgv4Addx7Zf7qHIOsCh3dIWl/Lvp9obdI9V0TZJY0qd60la/y3UsxXwQ+AA25sBU4Hzgff0kDdD2REREfG6QRW4LY+kYyTdWX66JnSdCGxaeq1OrOe3/QowBxhfyn9O0iWSrpR0n6R/lvQVSfMl/Y+ktUu+oyTdLen3kn5Rq3JP4DKqIGufbs37mKSbJf1B0idKPR2SNq21/5YStM0CZtteXNpp25fZ/l0t3/GSbgK+JGkjSbdJmkO1X2pEREQMUS0RuEnaFtgf2BaYAXxR0pZUQdDi0iM3q1uZkcA2wNW15C2ogq7tgG8DT9ieCswFDih5jgGm2N4K+FKt7H7AueVnv25NnAh8GPgUcLqkt1EFeHuXtkwAxtj+fWnDvD6+8mjbH7L9Q+Bk4CTb2wB/66NcREREtLGWCNyAHYCLbT9v+xmqnq/te8m7qaQFwGPAvbbvqp273vZztpcBzwJXlPQ7gEnl813ALyTtD3QCSBoPrA/cavtuYJikzWr1XmD7tdKLthSYDFwA7FXO71OO30TS2NJbuETSl2unzqt9nkEVBAL8vJfvHBEREUNAqwRuvb7PpAddc9w2Bj4sadfauZdqn1+rHb/GG69G+RjV/LptgQ5Jw6gCrzHAfZLupwri9q3V1f1leLb9Z+BZSZuX8l3B113A1iXTI6WtZwCjauWf61Z3XrYXERERLRO43QTsLmmkpFHATOBm4Bng7T0VsP0X4Gvlp19KkDbB9vXAV4B1gTWohkZ3sT3J9iSqoK4+XLqXKptQDZsuKennl+u/t2oVuQAACJNJREFUrfTUAXwH+EZ9/lu5Rm9upQy5Ug0XR0RExBDVEoGb7dup5pbNoQpkTrN9Rxny7JB0R/fFCcVFwDqSZvTzUqsBv5S0kGoe2reBscC7gI5ae5YAL0maVpLupQourwA+X1sBeyHwj9SGSW3PB/6lXGexpN9R9Q7Wh0frjgCOknQ7b+6Vi4iIiCEmW161gZFrjPL4z5yxyusdPudkFs27ZZXXGxEREcuVLa/a2fARwxnegO2pJo4bu8rrjIiIiJWXwK0NbLLxRnR0pGcsIiKi3bXEHLeIiIiISOAWERER0TKyOKENjF77HR6/4RarpK6J48ZyzZWXrJK6IiIiYqVkcUI763y5k85tDu87Yz8sbcAih4iIiFg1MlQaERER0SISuEVERES0iEEduEmaIOlXZRP2P0o6SdKIBl/z2fLnJEl31tK3lXRT2e3gHkk/kbS8rar6e73jJB29nPOzJS0sm9FfI2ncW71mREREtKZBG7hJEnAJcJntycAmVFs+Hf8W613heX2S1qPavuqrtjcF3gtcRS/7pK5i37W9ZdmM/krgGwNwzYiIiBiEBm3gBuwMvGj7pwC2XwWOAg6RNEfS68soJd0oaZqkNSWdWc7PlzSznD9Y0oWSrgCukTRK0nWS5pV9Tmf20ZbDgLNt/29pi21fZHuZpHUkXVZ6xW6VtGW55nGlLTdK+pOkI2rt/XrpubsW2LTnS1ZsP107XBPIMuCIiIghajCvKt0CmFtPsP20pAeoep72Br4p6d3AONtzJZ0AXG/7EElrA7eX4AhgBrCl7cdLr9vupb53ArdKuty9vxvlfcDZvZz7FjDf9m6SdgZ+Bkwp5zYDdqLqmVss6TRgS2BfYCrV/Z/X/Xt2J+l44EDgqVJfREREDEGDucdN9Ny7JOBGYK9yvDfVMCbAR4FZkhaUPKsD65dzv7H9eK2OEyQtBK4FxgPrrWQ7twd+DmD7emCMpLXKuV/bfsn2o8Aj5Ro7AJfafr70pl3e1wVsf932ROAc4Esr2c6IiIhocYM5cLsLmF5PkDQamAjMAR4rw5L7AOd1ZQH2sD2l/Kxve1E591ytqv2BdYFpZe7YMqogb3ltmdbLuZ5ektcVcL5US3uVN3o4V3a485fAHitZNiIiIlrcYA7crgPWkHQggKRhwPeBs2w/TxWsHQOsZfuOUuZq4PCysAFJU3upey3gEdudknYC3tNHW04BDpL0ga4ESQdIehdwE1UgiKQdgUe7zUvr7iZgd0kjJb0d+NTyLixpcu3w08A9fbQ1IiIi2tSgDdzKfLPdgb0kLQH+ALwIHFuyXEQ1V+yCWrHZwHBgYXmVx+xeqj8HmC6pgyroWm4wZHtZudb3yqKCRVRDnk8Dx5W6FgInAgf1Udc84HxgAXAxcPPy8gMnSrqz1P9R4Mg+8kdERESbyl6lbSB7lUZERLSV7FXazjbZeCM6Om5pdjMiIiKiwRK4DRKSTgU+2C35pK732EVEREQkcBskbB/W7DZERETE4JY5bm1A0jPA4ma3Ywh7J/BosxsxhOX+N0/ufXPl/jdPo+/9o7Y/3tOJ9Li1h8W2p/edLRpBUkfuf/Pk/jdP7n1z5f43TzPv/aB9HUhEREREvFkCt4iIiIgWkcCtPZze7AYMcbn/zZX73zy5982V+988Tbv3WZwQERER0SLS4xYRERHRIhK4tThJHy/7p94raVaz29POJE2UdIOkRZLuknRkSV9H0m8kLSl/vqPZbW1nkoZJmi/pynK8gaTbyv0/X9KIZrexXUlaW9JFku4pvwcz8vwPDElHlX937pR0rqTV8+w3jqQzJT1S9j3vSuvxWVflR+X/4YWStm5k2xK4tTBJw4BTgU8AmwP7Sdq8ua1qa68A/2r7vcB2wGHlfs8CrrM9GbiuHEfjHAksqh1/G/hBuf9PAP/UlFYNDScBV9neDNiK6u8hz3+DSRoPHAFMt/0+YBiwL3n2G+ksoPt71Hp71j8BTC4/nwdOa2TDEri1tm2Be23/yfbLwHnAzCa3qW3Zftj2vPL5Gar/tMZT3fOzS7azgd2a08L2J2kC8A/AT8qxgJ2Bi0qW3P8GkTQa+BBwBoDtl20/SZ7/gbIaMFLSasAawMPk2W8Y2zcBj3dL7u1Znwn8zJVbgbUlvbtRbUvg1trGA0trxw+WtGgwSZOAqcBtwHq2H4YquAPGNq9lbe+HwDHAa+V4DPCk7VfKcX4HGmdD4G/AT8tQ9U8krUme/4az/RDwPeABqoDtKWAuefYHWm/P+oD+X5zArbWph7QsE24wSaOAi4Ev23662e0ZKiR9EnjE9tx6cg9Z8zvQGKsBWwOn2Z4KPEeGRQdEmUs1E9gAGAesSTU8112e/eYY0H+HEri1tgeBibXjCcBfmtSWIUHScKqg7Rzbl5TkZV3d4uXPR5rVvjb3QeDTku6nmhawM1UP3Npl+AjyO9BIDwIP2r6tHF9EFcjl+W+8XYD7bP/NdidwCfD35NkfaL096wP6f3ECt9Y2B5hcVhaNoJqsenmT29S2ynyqM4BFtv+jdupy4KDy+SDgVwPdtqHA9tdsT7A9iepZv972/sANwJ4lW+5/g9j+K7BU0qYl6SPA3eT5HwgPANtJWqP8O9R17/PsD6zenvXLgQPL6tLtgKe6hlQbIS/gbXGSdqXqdRgGnGn7+CY3qW1J2h64GbiDN+ZYHUs1z+0CYH2qf2D3st19UmusQpJ2BI62/UlJG1L1wK0DzAcOsP1SM9vXriRNoVoYMgL4E/BZqg6APP8NJulbwD5Uq9vnA5+jmkeVZ78BJJ0L7Ai8E1gGfBO4jB6e9RJMn0K1CvV54LO2OxrWtgRuEREREa0hQ6URERERLSKBW0RERESLSOAWERER0SISuEVERES0iARuERERES0igVtEREREi0jgFhEREdEiErhFREREtIj/A5QU9AcPhsUgAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 648x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from statlearning import plot_feature_importance\n",
    "\n",
    "plot_feature_importance(gb.best_estimator_, predictors)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### XGBoost\n",
    "\n",
    "[XGBoost](https://xgboost.readthedocs.io/en/latest/) is a state-of-art gradient boosting library that is very popular among [Kaggle](https://www.kaggle.com/) users. The easiest way to get started with XGBoost is to use the [Scikit-Learn API](https://xgboost.readthedocs.io/en/latest/python/python_api.html#module-xgboost.sklearn) provided by the package. That makes the syntax identical to what we did above, except that we call the [XGBRegressor](http://xgboost.readthedocs.io/en/latest/python/python_api.html#xgboost.XGBRegressor) class from the XGBoost package.  (Note that the <TT>xgboost</TT> package needs to be installed, see the beginning of the tutorial.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Best parameters found by randomised search: {'subsample': 0.6, 'n_estimators': 1000, 'max_depth': 2, 'learning_rate': 0.1} \n",
      "\n",
      "CPU times: user 16.8 s, sys: 127 ms, total: 16.9 s\n",
      "Wall time: 11.6 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "model = xgb.XGBRegressor()\n",
    "\n",
    "tuning_parameters = {\n",
    "    'learning_rate': [0.01, 0.05, 0.1],\n",
    "    'n_estimators' : [250, 500, 750, 1000, 1500],\n",
    "    'max_depth' : [2, 3, 4],\n",
    "    'subsample' : [0.6, 0.8, 1.0],\n",
    "}\n",
    "\n",
    "gb_search = RandomizedSearchCV(model, tuning_parameters, n_iter = 1, cv = 5, return_train_score=False, n_jobs=4,\n",
    "                              random_state = 20)\n",
    "gb_search.fit(X_train, y_train)\n",
    "\n",
    "xbst = gb_search.best_estimator_\n",
    "\n",
    "\n",
    "print('Best parameters found by randomised search:', gb_search.best_params_, '\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnoAAAF1CAYAAAB/FEdEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdebid093/8fenkUgiEvOQCDEkZg5JEDNFRyWmVBWpqvJoPLSqHp1SnqKlAwlVU0OLKmJo9EdMITHlnEQkIQklIcVjKEUGyRHf3x/32nLbzj5DcmJn7/15Xde5cu91r+G7d/75Xmvd616KCMzMzMys+nyu3AGYmZmZ2YrhRM/MzMysSjnRMzMzM6tSTvTMzMzMqpQTPTMzM7Mq5UTPzMzMrEo50TMzawNJV0v6fSvrbiEpJG3QTJ0Jks5ejnjGSvr+srY3s+q2SrkDMDNrb5LuBN6NiOOauPcQ8ExEfG9Z+o6IE5c3vvYUEQeVO4amSPoLMC8iTi53LGa1zDN6ZlaN/ggcIWmNfKGkvsA+wJVt7VBSB0lqp/iqln8ns5WLEz0zq0b3AG8CxxaVnwQ8ERFTAST9StJsSfMk/VPSsELF3LLrtyTNBBYAa0v6i6QrcvVK9pHzFUnPS/qPpNslrVsqcEl9JI2W9H+SXpX0B0mrNVP/46XfXMzHSZohab6kMZLWkPRrSW9Kek3Sd3PtT5Q0U9I5aczXU91VcnXqJD0k6R1JL6S6HZr5nX4KDAG+nX6XeanuTpIekfRvSW9LulvSprlx/iLpT5KuSb/VvyR9YgZV0n6SHk2xvCnp6ty9HSTdJ+ktSS9L+qWkjqV+O7Na4ETPzKpORHwEXA18p1AmqRNwPJ+czZsO7A6sDpwMXCTp80XdHQ3sDXQH3mliuNb08U1gD2AToAMwqqm4JXUFHgKeBjYFtkv//q6579uEwcCgNN4WwERgBrAB2W8yQlKvXP3NgfXTWHsChwFnpJjWBO5Lf+sDB5MlzKcVjZn/nX4J3AxcExHdIqJbqhNkSeCGwGbAIuD6on6GAKOBtYDvA5dL2ijFshPw/4Ar0nfZBPhLurcB8HAad0Oy/5MvAT9s7Y9mVo2c6JlZtboG2FrSrunzYKAj8LdChYj4c0S8Fpn7yWYCi5O04RHxRkQsioglxYO0sY93gbOAL0tar4mYvwY0RsQvImJhRLwN/Bw4to3LoedGxH8i4i3gH8DCiPhTRCyJiDHAPKAuV78RODuN+TxwMfCtdO9gYD5wQUQsjohngYuA4mcVm/2dACJiSkQ8nPr5D3AusIekzrlq90XE3RHxUUT8LY29Y7p3CnB7+s0XRcSCiBiX7g0F6iPi6ohojIh/Ab8CPvWcplkt8WYMM6tKEfGqpLvJZp+eTP/+OSIWFOpIOh34NtA7FXUlW/LNm9PcOMvQR+F6I+CNonqbAptK+k9+iPS3HvB6c7HkvJa7XlD0uVC2eu7z6xGxsCjGjdJ1b2B2RETu/gss/b75Ns1Kz0j+GtglN76AtYFXmogdskSvULcP8HiJ7jcF9in67T4HfNRSXGbVzDN6ZlbN/ggMSUt++5FbtpW0D9kS43eAtSNiDbLZr+KZs5KJQhv66NPE9b+a6PIl4NmIWCP31yMiOkdEa5O8ZbF+0axan1x8c4E+RTOKm6XyvOLfqanf7Uqy5e/tI6I72VIvfPr3KmUO0LfEvZeAe4p+u+7p/8SsZjnRM7Nqdi/wFnAb8HhETM/d6w4sIZt9C0kHA219VUlr+/iZpHUl9QAuJEtIimfzAO4Cukn6kaRuymwk6dA2xtVWHYELJHWWtAXwA+C6dO/vZDNqZ0nqJGlrsufermmhz/8DNi9KELuTLRu/mzak/KKNcV4BHCbpGymWrpL2TfdGAYMkHZ++x+ckbS7pC20cw6yqONEzs6qVNmVcRbasV/xKlX8ANwGTyBK1Q4E72zhEa/u4EXgMeJls9mpoiXjnkc087gjMAt4l2wSxfRvjaqsXyOKfk+K8C/htiukdsuT1y2RLzf8gS/IuaaHPK4E1gLdzy6mnA/sD7wHjyJLIVouIyWTPDJ6W4n0J+Ea692rq+4hU/jZZgt+nLWOYVRt98rELMzOrJen1JWdGxFbljsXM2p9n9MzMzMyqlBM9MzMzsyrlpVszMzOzKuUZPTMzM7Mq5UTPzMzMrEr5ZIwa9MUvfjHuueeecodhZmZm7aPkS8c9o1eD3nrrrXKHYGZmZp8BJ3pmZmZmVcqJnpmZmVmV8utValD3NdaMXpttW+4wzMzMakLvnusxdszoFTlEyWf0vBmjBjUubqRx4LByh2FmZlYT5taPKNvYNb90K2l9STdKelHSJEmPSxpcxni+JKlB0gxJMyVdXK5YzMzMrLLVdKInScAdwCMRsVlE9Ae+DmzUyvYd2jme7YCRwDcjYmtgO+DFNrT3DK2ZmZl9rKYTPWB/YHFEXFEoiIiXImKEpD6SxkuanP52B5C0r6SHJN0ITEtld6TZwGcknVToS9K3JT0naZykqySNTOXrSrpNUn362yM1OQv4ZUTMTLF8GBGXpzYHS3pS0lOS7pe0fiofLulKSWOB6yVtK2mipCmSpkrqu8J/RTMzM1sp1foM0LbA5BL33gAOjIgPUrJ0EzAg3dsF2C4iZqfPJ0TE25K6APWSbgNWBX4K7Ay8DzwIPJ3qXwL8LiImSNoYuBcozOD9pkQ8E4DdIiIknUiWFP4g3esP7BkRCyWNAC6JiBskdQLaddbRzMzMKketJ3qfIOkyYE9gMXAAMFJSHbAE6JerOjGX5AGclnuurzfQF9gAeDgi3k5935Lr4wBgm2zlGIDuklZvIbyNgJslbQh0AvLj3xURC9P148CPJW0EjI6I51vx1c3MzKwK1frS7TNkM24ARMSpwOeBdYEzgNeBHclm8jrl2s0vXEjalyxxGxQROwJPAZ1pZqsz2e8+KCLq0l+viHg/xdO/RJsRwMiI2B74bhrjU/FExI3A14CFwL2S9m8mDjMzM6titZ7oPQh0lnRKrqxr+rcH8FpEfAQcS+kl0B7AOxGxQNJWwG6pfCKwj6Q10yaJw3NtxgLfK3xIs4YAFwHnSOqXyj8n6fu5cV5J18eX+kKSNgNejIhLgbuAHUrVNTMzs+pW04leZG+LPpQsIZstaSJwHfAj4HLgeElPkC25zi/RzT3AKpKmAucBT6S+XwHOB54E7geeBd5NbU4DBqTNEs8CJ6c2U4HTgZskzQCmAxumNsOBWySNB5o7rHYIMF3SFGAr4PrW/yJmZmZWTXwyxgokqVtEzEszercD10bE7eWOq0vXbtHr2GvKHYaZmVlN6Fg/ghmTJ6zIIXwyRpkMl3QA2fN0Y8ne2Vd2HTt1pGMZ39JtZmZWS3r3XK9sY3tGrwYNGDAgGhoayh2GmZmZtY+SM3o1/YyemZmZWTVzomdmZmZWpZzomZmZmVUpP6NXg7qvsWb02mzbcodh1m5691yPsWNGlzsMM7Ny8a5bW6pxcSONA4eVOwyzdjPXu8jNzJrkpdsWSFoiaUrur4+kAZIubccx5khap736MzMzMwPP6LXGwoioKyqbA3zq/SSSVomIDz+TqMzMzMxa4Bm9ZSBpX0lj0vVwSVdKGgtcL6mDpIsk1acjzr6ba/OIpNslPSvpCkmf+v0l3SFpkqRnJJ2UK/+ipMmSnpb0QCpbTdK1aaynJB2SyreVNDHNQE6V1Pcz+WHMzMxspeIZvZZ1SefGAsyOiMFN1OkP7BkRC1Ny9m5EDJS0KvBoSgIBdgG2AV4iOyP3MODWor5OiIi3JXUB6iXdRpaQXwXsHRGzJa2V6v4YeDAiTpC0BjBR0v1kZ+deEhE3SOoEdGiPH8LMzMwqixO9ljW1dFvsrohYmK4PAnaQdET63APoCywGJkbEiwCSbgL25NOJ3mmSCslk79R2XeCRiJgNEBFv58b6mqQz0+fOwMbA48CPJW0EjI6I59v0jc3MzKwqONFrH/Nz1wKGRcS9+QqS9gWK32UTTdQ5ABgUEQskjSNL3tRE28JYh0fErKLyGZKeBL4C3CvpxIh4sE3fyMzMzCqen9Frf/cCp0jqCCCpn6TV0r1dJG2ans0bAkwoatsDeCcleVsBu6Xyx4F9JG2a+iws3d4LDJOkVL5T+ncz4MWIuBS4C9hhRXxRMzMzW7k50Wt/VwPPApMlTQf+yNKZ08eBC4HpwGzg9qK29wCrSJoKnAc8ARARbwInAaMlPQ3cnOqfB3QEpqaxzkvlQ4Dp6dnCrYDr2/tLmpmZ2crPJ2N8RtKy7JkR8dVyx9Kla7fodew15Q7DrN10rB/BjMnFE+RmZjXDJ2PYUh07daSjTxKwKtK753rlDsHMbKXkGb0aNGDAgGho+NT7ns3MzKwylZzR8zN6ZmZmZlXKiZ6ZmZlZlXKiZ2ZmZlal/IxeDeq+xprRa7Ntyx2G2XLp3XM9xo4ZXe4wzMxWBt51a0s1Lm6kceCwcodhtlzmeue4mVmLvHRrZmZmVqUqLtGTNC93/WVJz0vaWNLJko5L5UMl9Wyhn6GSRrZjXIdKmipppqTpko5Yjr76pJMumrr3BUlT0t88SbPS9fW5OpdIeiUdtWZmZmY1qmKXbiV9HhgBHBQRLwNX5G4PJTtm7NXPKJYdgYuBAyNidjqT9n5JsyNiUnuOFRH3kp1xi6RxZKdtfPxSvJTcDQbmAnsD49pzfDMzM6scFTnjI2kv4CrgKxHxQiobLunMNJM2ALghzXR1kTRQ0mOSnpY0UdLqqaueku5Js4K/zvV/kKTHJU2WdIukbql8jqRfpPJpkrZKTc4Ezo+I2QDp3/OBH6R24yQNSNfrSJqTrvtIGp/6myxp93b4efYjS3L/ABzdDv2ZmZlZharERG9V4E7g0IiYWXwzIm4FGoBjIqIOWALcDPx3ROwIHAAsTNXrgCHA9sAQSb0lrQP8BDggInZOfX0/N8RbqfwPZAkewLZA8cxdA7BNC9/lDbJZwJ1THJe29OVb4WjgJuB24KuSOrZDn2ZmZlaBKjHRawQeA77dyvpbAq9FRD1ARLwXER+mew9ExLsR8QHwLLAJsBtZgvaopCnA8am8oPA+h0lAn3QtoPg9NSW3Oud0BK6SNA24hZYTw2ZJ6gR8GbgjIt4DngQOWp4+zczMrHJV4jN6HwFHkT0Dd05EnN9C/aaSsIJFueslZL+HgPsiotSy56Ki+gDPkC0XT83VK8wGAnzI0qS6c67OGcDrwI7p/gfNfZFW+CLQA5gmCaArsAC4ezn7NTMzswpUiTN6RMQC4KvAMZKamtl7Hyg8hzeT7Fm8gQCSVpfUXIL7BLCHpC1S/a6S+rUQ0sXA/0jqk9r0AU4HLkr35wD903V+N24PstnGj4BjgQ4tjNOSo4ETI6JPRPQBNgUOktR1Ofs1MzOzClSRiR5ARLxNNoP1E0mHFN0eBVyRll47kD3/NkLS08B9fHJWrbjfN8l27d4kaSpZ4rdVqfqpzRTgR8DfJT0HPAecEhGzUpWLgVMkPQask2t6OXC8pCeAfsD8lr53KSmZ+wK52buImA9MAA5e1n7NzMyscvkItBVA0oXArsAXImJxueMp5iPQrBr4CDQzs4+V3BfgRK8GDRgwIBoaGlquaGZmZpXAZ91WIklfAH5VVDw7IgaXIx4zMzOrLE70VmL5UzDMzMzM2qpiN2OYmZmZWfP8jF4N8mYMq0TefGFmVpKf0bOlGhc30jhwWLnDMGuTufUjyh2CmVnF8dKtmZmZWZWqqERP0vqSbpT0oqRJkh6XVLYdqJK+JKlB0gxJMyVd3E79jpJ0RIl7t0uaIumfkt5N11Mk7Z7uryupUdJ32yMWMzMzq1wVk+gpO7z1DuCRiNgsIvoDXwc2amX75T1erLi/7YCRwDcjYmtgO+DF9hyjKRExOCLqgBOB8RFRl/4eS1WOJDvNo9RZvWZmZlYjKibRA/YHFkfEFYWCiHgpIkZI6iNpvKTJ6a8wu7WvpIck3QhMS2V3pNnAZySdVOhL0rclPSdpnKSrJI1M5etKuk1SffrbIzU5C/hlRMxMsXwYEZenNptIekDS1PTvxql8lKRLJT2WZiWPSOWSNFLSs5LuBtZbjt/paOAHwEaSei1HP2ZmZlbhKinR2xaYXOLeG8CBEbEz2bm2l+bu7QL8OCK2SZ9PSLOBA4DTJK0tqSfwU2A34EA+ebbtJcDvImIgcDhwdSrfDphUIp6RwPURsQNwQ1E8GwJ7Al8FLkxlg4Etge2B7wC7l+i3WZJ6AxtExETgb2S/hZmZmdWoit11K+kysoRpMXAAMFJSHbAE6JerOjEiZuc+n5Z7rq830BfYAHg4It5Ofd+S6+MAYJts5RiA7pJWbyG8QcBh6frPwK9z9+6IiI+AZyWtn8r2Bm6KiCXAq5IebKH/Ur5OluAB/BW4BvjtMvZlZmZmFa6SEr1nyGbUAIiIUyWtAzQAZwCvAzuSzVJ+kGs3v3AhaV+yxG1QRCyQNA7oTDPvn0n9DYqIhflCSc8A/YGnWxF7/mWFi/LdlKizrI4G1pd0TPrcU1LfiHi+Hfo2MzOzClNJS7cPAp0lnZIr65r+7QG8lmbKjgVKbbzoAbyTkrytyJZqASYC+0haU9Iq5BJKYCzwvcKHNGsIcBFwjqR+qfxzkr6f7j1GNrsGcAwwoYXv9gjwdUkdJG0I7NdC/U+RtCWwWkT0iog+EdEHuCAXh5mZmdWYikn0IjvC41CyhGy2pInAdcCPgMuB4yU9QbbkOr9EN/cAq0iaCpxHtjuViHgFOB94ErgfeBZ4N7U5DRiQNlY8C5yc2kwFTgdukjQDmE72/F2hzbfSOMcC/93C17sdeJ5sw8gfgIdb9aN80tGpn7zb8O5bMzOzmuUj0BJJ3SJiXprRux24NiKKE6eq4CPQrBL5CDQzs5J8BForDJd0ANkze2PJ3tlXlfptsTkNDS2tJpuZmVmlc6KXRMSZ5Y6hmKTbgU2Lin8UEfeWIx4zMzOrLE70VmIRUbbj3czMzKzyVcxmDDMzMzNrG2/GqEHejGGfFW+gMDP7THgzhi3VuLiRxoHDyh2G1YC59SPKHYKZWU3z0m07kLRE0pTcX58W6s9Jp3ogaV76t4+khan905IeSy9Bbq6fPpK+kfs8VNLI5f9GZmZmVg2c6LWPhRFRl/ubs4z9vJDa70j2MuhzWqjfB/hGC3XMzMysRjnRW0GKZ9ckjUln7bZWd+Cd1LaPpPGSJqe/3VOdC4G90izgGamsp6R7JD0v6dft8V3MzMysMvkZvfbRRdKUdD17OV6LsnnqZ3Wyc3x3TeVvAAdGxAeS+gI3AQOAs4EzI+KrkCWXQB2wE7AImCVpRETMXcZ4zMzMrII50WsfCyOirh36eaHQj6QhwJXAF4GOwEhJdcASsvN8S3kgIt5NfTwLbAI40TMzM6tBTvRWnA/55NJ45za2vwv4U7o+A3gd2DH1+UEz7Rblrpfg/2MzM7Oa5Wf0Vpw5QJ2kz0nqDezSxvZ7Ai+k6x7AaxHxEXAs0CGVv0+2zGtmZmb2KZ7tWXEeBWYD04DpwORWtCk8oydgMXBiKr8cuE3SkcBDwPxUPhX4UNLTwCjS5g0zMzMz8MkYNalL127R69hryh2G1YCO9SOYMXlCucMwM6t2JU/G8NKtmZmZWZXy0m0N6tipIx19NJV9Bnr3XK/cIZiZ1TQnejWo3xab09Dg5TQzM7Nq56VbMzMzsyrlRM/MzMysSnnXbQ3qvsaa0Wuzbcsdhq3kevdcj7FjRpc7DDMza1nJXbd+Rq8GNS5upHHgsHKHYSu5ud6wY2ZW8bx0a2ZmZlalnOiZmZmZVamqTvQkLZE0Jfd3dgv1z1nGcTpJ+r2kFyT9U9IYSRsvW9QgabikM0vcuyx9l2clLcx9tyPS/VUkvSXpgmUd38zMzKpDtT+jtzAi6tpQ/xzg/LYMIKlDarM60C8ilkj6FnCnpP4R8VFb+mtJRJyaxu0DjGni+x0EzAKOknROeLeNmZlZzarqGb2mSOohaZakLdPnmyR9R9KFQJc0O3ZDuvdNSRNT2R9TUoekeZLOlfQksAfwLeCMiFgCEBF/AuYBB0jqI2l6bvwzJQ1P19+RVC/paUm3SeraDl/xaOAS4GVgt3boz8zMzCpUtSd6hcSt8DckIt4FvgeMkvR1YM2IuCoizibNAEbEMZK2BoYAe6RZsyXAManf1YDpEbEr8B/g5Yh4r2jsBmCbFuIbHREDI2JHYAbw7eX5spK6AJ8HxgA3kSV9ZmZmVqNqcuk2Iu6TdCRwGbBjibafB/oD9ZIAugBvpHtLgNvStYCmlkdLvtMmZztJ/wusAXQD7m1Fm+Z8FXgoIhZIug34qaSPZxrNzMystlR7otckSZ8DtgYWAmsB/2qqGnBdRPxPE/c+yCVP/wQ2kbR6RLyfq7MzcCvwIZ+cOe2cux4FHBoRT0saCuzb9m/zCUcDe0iakz6vDewH3L+c/ZqZmVkFqval21LOIFsqPRq4VlLHVN6Yu34AOELSegCS1pK0SXFHETEfuA74be4ZvuOAD4BHgdeB9SStLWlVslm3gtWB19KYx7AcJHUH9gQ2jog+EdEHOBUv35qZmdWsap/R6yJpSu7zPcC1wInALhHxvqRHgJ8APweuBKZKmpye0/sJMDbNADaSJU4vNTHO/wAXAbPSc3JvAoPSjtdGSecCTwKzgZm5dj9N5S8B08gSv2V1GPBgRCzKld0J/FrSqkXlZmZmVgN81m07k7QBWUJ5eURcWe54muKzbq01fNatmVnFKLkvwIleDRowYEA0NDSUOwwzMzNrHyUTvWpfuq1oki4je09f3iXpPX1mZmZmzXKitxIrnIJhZmZmtixqddetmZmZWdXzM3o1yJsxysubHMzMrJ35GT1bqnFxI40Dh5U7jJo1t35EuUMwM7Ma4aVbMzMzsypVlYmepI0k3SnpeUkvSLpEUqcVPOa89G8fSdNz5XtKmihppqRZkpZrg0VhnBbqdJf0iqSRyzOWmZmZVbaqS/QkCRgN3BERfYF+QDfgl8vZb5uXudPLk28ETo6IrchelXKCpMHLE0srnAc8vILHMDMzs5Vc1SV6wP7AB4V3zUXEErKzbU+QVC/p410IksZJ6i9pNUnXpvtPSTok3R8q6RZJfyc7Cq2bpAckTZY0rVCvGacCoyJicorlLeAs4Iep/1GSjsjFU5gVbOs4H5PUH1gfGNvaNmZmZladqjHR2xaYlC+IiPeAl4ExwFEAkjYEekbEJODHZOfEDgT2Ay6StFpqPgg4PiL2Bz4ABkfEzqneb9IMYqtjARqAbVr4Dm0dh/SdPgf8hpRImpmZWW2rxkRPQFPvjBEwDjgyfT4KuCVdHwScLWlKqtMZ2Djduy8i3s71cb6kqcD9QC+y2bO2xtKa79CWcQr+C/hHRMxdhjHNzMysylTj61WeAQ7PF0jqDvQG6oF/S9oBGAJ8t1AFODwiZhW12xWYnys6BlgX6B8RjZLmkCWFzcUyALgrV9afbFYP4ENSsp1m7AobRto6TsEgYC9J/0X2XGInSfMi4uxWtDUzM7MqU40zeg8AXSUdByCpA9ly5qiIWAD8lew5uR4RMS21uRcYVlgelbRTib57AG+k5Gs/YJMWYrkMGCqpLvW7NtmmkPPS/TlkiR/AIUDHZRwHgIg4JiI2jog+wJnA9U7yzMzMalfVJXqRHfUxGDhS0vPAc2TPvJ2TqtwKfB34W67ZeWRJ1tT0apTzaNoNwABJDWSzbjNbiOU14JvAlZJmAa8Cl0ZEYUfsVcA+kiYC+dnDNo1jZmZm1hQfgfYZSu/QOxnYOyLeKVccPgKtvHwEmpmZtbOSGzad6NWgAQMGRENDQ8sVzczMrBL4rNtqIml74M9FxYsiYtdyxGNmZmYrJyd6FShtIqkrdxxmZma2cqu6zRhmZmZmlvEzejXImzFWDG+yMDOzMvEzerZU4+JGGgcOK3cYVWdu/Yhyh2BmZvYJXro1MzMzq1JVk+hJ2kjSnZKel/SCpEskdWq55XKNOS/92ye9aLlQvoukRyTNkjRT0tWSurbDeMMlndnM/fMkTZU0RdJYST2Xd0wzMzOrXFWR6KWjy0YDd0REX6Af2Vmvv1zOftu8tC1pfeAW4EcRsSWwNXAPsPryxNJKF0XEDhFRB4wBfvYZjGlmZmYrqapI9ID9gQ8i4k8AEbEEOAM4QVK9pI93HkgaJ6m/pNUkXZvuPyXpkHR/qKRbJP0dGCupm6QHJE2WNK1QrxmnAtdFxOMploiIWyPidUlrSbojzbo9IWmHNObwFMs4SS9KOi0X74/TzOD9wJbNDRwR7+U+rgZ4p42ZmVkNq5bNGNsCk/IFEfGepJfJZraOAn4uaUOgZ0RMknQ+8GBEnCBpDWBiSqYABgE7RMTbaVZvcOpvHeAJSXdF6e3K2wHXlbj3C+CpiDhU0v7A9Sx9H95WwH5kM3+zJP0B2IHsXN6dyP6vJhd/z2KSfgkcB7yb+jMzM7MaVS0zeqLp2SsB44Aj0+ejyJZVAQ4CzpY0JdXpDGyc7t0XEW/n+jhf0lTgfqAXsP4yxrkn6USLiHgQWFtSj3Tv7ohYFBFvAW+kMfYCbo+IBWm27q6WBoiIH0dEb+AG4HvLGKeZmZlVgWpJ9J4BBuQLJHUHegP1wL/TMukQ4K+FKsDhEVGX/jaOiBnp3vxcV8cA6wL907Nvr5Mlhc3F0r/Evabec1NIUBflypawdLZ1WZdfbwQOX8a2ZmZmVgWqJdF7AOgq6TgASR2A3wCjImIBWXJ3FtAjHR8GcC8wLG3kQNJOJfruAbwREY2S9gM2aSGWkcDxkj4+d1bSNyVtADxCljgiaV/graLn6oo9AgyW1EXS6sDBzQ0sqW/u49eAmS3EamZmZlWsKhK99LzcYOBISc8DzwEfAOekKreSPev2t1yz84COwNT0apTzSnR/AzBAUgNZktZs8hQRr6exLk6bKGaQLcG+BwxPfU0FLgSOb6GvycDNwBTgNmB8c/WBCyVNT/0fBPx3C/XNzMysivkItBrkI9BWDB+BZmZmZeIj0GypfltsTkPDhHKHYWZmZiuYE70KJOkyYI+i4ksK7xE0MzMzAyd6FSkiTi13DGZmZrbyq5mXJC0AACAASURBVIrNGGZmZmb2ad6MUYO8GaPtvNHCzMxWYt6MYUs1Lm6kceCwcodRUebWjyh3CGZmZm1WU0u3ktaWNCX9/Z+kV3KfOzVRfy1JJ7ei31Uk/SddbyFpYerzaUmPFr3IeFlj31/SbrnPW0t6OI0zI52Ni6QDJL2b+173Lu/YZmZmVplqakYvIv4N1AFIGg7Mi4iLm2myFnAycEUbh5qVjktD0qnA2cC32xzwJ+0PvAU8kT6PBH4dEXen0z22y9V9KCIOXc7xzMzMrMLV1IxecySdlU6VmC6psK55IbBlmhm7UFJ3SQ9KmixpqqSvtqLr7sA7aYztJdWn/qZK2izNAE6XdK2kZyRdL+kLkh6T9JykAZI2B04Efpja7g5sCPwLspNBcke7mZmZmQE1NqNXiqRdyI432wXoAEyU9DDZTNwWudm5jsAhEfG+pPWAR4ExTXS5paQpZEneqkDh3Nv/Ai6OiJslrUr28ORGwJbAUWTHq00GFkXE7pIOB86OiCMkXU12Nu7vUyy/BR6R9CgwFvhTRLybxtkvjQ/w14i4sF1+KDMzM6sontHL7AXcFhELIuJ94A5gzybqCfhVOkt2LNBb0jpN1JsVEXURsRlwFkuXfh8DfiLpLKB3RHyQyv8ZEc9GxEfAs8D9qXwa0KepgCPiamAbsnN8Pw88nnvO8KE0fp2TPDMzs9rlRC9TcltykeOAHsDOaZbvLaBzC23uAvYGiIg/A4OBRcB9kvZOdRbl6n+U+/wRzcy6RsQrEXFtRBxM9n+5dSu/h5mZmdUAJ3qZR4DBkrpI6gYcAowH3gdWz9XrAbwRER9KOhDo1Yq+9wReAJC0WUT8MyIuAe4GdmhDjJ+IRdIXJa2SrnsCawKvtqE/MzMzq3J+Rg+IiImSbgLqU9EfCpsbJDVImkaWmP0W+LukBrJn6Z4v0WXhGT2Rzc6dlMq/IelooJEsKfsJ0NTSb1PuBG6RdBhwKvAl4BJJHwABnB4Rb2YbcM3MzMx8MkZN6tK1W/Q69ppyh1FROtaPYMbkCeUOw8zMrCk+GcOW6tipIx190kOb9O65XrlDMDMzazMnejWo3xab09Dg2SkzM7Nq580YZmZmZlXKiZ6ZmZlZlXKiZ2ZmZlalvOu2BnVfY83otdm25Q5jpde753qMHTO63GGYmZm1xLtubanGxY00DhxW7jBWenO9M9nMzCqcl27NzMzMqlTZEz1JSyRNkfS0pMmSdm+HPuskfTn3eaikN9M4hb9tlneccpHUV9IYSS9ImiTpody5ucV150hq7ekbZmZmVkVWhqXbhRFRByDpC8AFwD7L2WcdMAD4R67s5oj43nL22+4krRIRH7ahfmey49jOjIi7Utl2ZN/3kRUTpZmZmVWiss/oFekOvAMgaUNJj6TZt+mS9krl8yT9Ks1k3S9pF0njJL0o6WuSOgHnAkNS2yGlBpM0OPWhNN5zkjZIM4B3SrpH0ixJP8+1+X6KZ7qk01PZapLuTrOS0wtj5mfTJA2QNC5dD5d0paSxwPWSOki6SFK9pKmSvtvMb3QM8HghyQOIiOkRMSr1vbaksZKekvRHmnlA08zMzKrbyjCj10XSFKAzsCGwfyr/BnBvRPxSUgegaypfDRgXET+SdDvwv8CBwDbAdRFxl6SfAQMKM3iShpIlfnvmxh0UEbdLOhw4Ffgi8POI+D9JALsA2wELgHpJdwMBfAvYlSyBelLSw8BmwKsR8ZU0Xo9WfO/+wJ4RsVDSScC7ETFQ0qrAo5LGRsTsJtptC0xupt+fAxMi4lxJXwFOakUsZmZmVoVWhkQvv3Q7iGyGazugHrhWUkfgjoiYkuovBu5J19OARRHRKGka0KeZcUot3Q4DpgNPRMRNufL7IuLfKa7RwJ5kid7tETE/V75XiudiSb8CxkTE+FZ877siYmG6PgjYQdIR6XMPoC/QVKL3CSnZ7Qs8FxGHAXsDhwFExN2S3mlFLGZmZlaFVqql24h4HFgHWDciHiFLWl4B/izpuFStMZa+/O8jYFFq+xHLlrj2Sv2sLyn/exS/YDAosQwaEc+RzdBNAy5IM4oAH7L0N+5c1Gx+7lrAsIioS3+bRsTYEvE+A+ycG3swMBRYq5nYzczMrAatVImepK2ADsC/JW0CvBERVwHXkEtuWuF9YPVWjLcK8CeyZeIZwPdztw+UtJakLsChwKNkmx0OldRV0mrAYGC8pJ7Agoj4C3BxLtY5ZAkgwOHNhHIvcEqavURSv9R/U24E9pD0tVxZ19z1I2TP8SHpS8CazYxrZmZmVWxlWLotPKMH2czW8RGxRNK+wA8lNQLzgONKddCEh4CzU78XpLLiZ/T+CzgAGB8R41PdwrN4ABOAPwNbADdGRAOApFHAxFTn6oh4Ku0WvkjSR0AjcEq6/wvgGknnAE82E+/VZMvOk5U9IPgmWXL5KemZvq8Cv5X0e+B1ssT2f3Nj3iRpMvAw8HIz45qZmVkV8xFoTUibNwaUeKav4vkItNbxEWhmZlYhfASaLdVvi81paJhQ7jDMzMxsBXOi14T0TrpR5YxB0vZkS8d5iyJi13LEY2ZmZpXHid5KKiKmkZ3wYWZmZrZMVqpdt2ZmZmbWfrwZowZ5M0bTvPnCzMwqlDdj2FKNixtpHDis3GGsdObWjyh3CGZmZu3KS7dmZmZmVapqEz1JSyRNkTRd0i2SurbcapnHGippZLoeLumVNPbzkkZL2qaZtudKOqCd4jg/nbdb+LyJpBclrdEe/ZuZmVllqdpED1iYzo3dDlgMnPwZjv27NHZf4GbgQUnrFleS1CEifhYR97fTuOcBh0jaOn2+BPhpRPynnfo3MzOzClLNiV7eeLKjzJB0h6RJkp6RdFIqO0XSrwuV0wzdiHT9TUkT0wzdHyV1SOXfkvScpIeBPUoNHBE3A2PJztNF0hxJP5M0AThS0ihJR0j6kqS/5WLYV9Lf0/VBkh6XNDnNTnYrMdZCsvN6L0/n3K4eETcs+89mZmZmlazqEz1JqwBfAqalohMioj8wADhN0trArcBhuWZDgJvTzNgQYI+IqAOWAMdI2pDsTNk9gAOBkkuzyWRgq9znDyJiz4j4a67sPmA3SasVxbAO8BPggIjYGWggS+aaFBH/AN4Gric7z9fMzMxqVDXvuu0iaUq6Hg9ck65PkzQ4XfcG+kbEE+lZtt2A54EtgUeBU4H+QL0kgC7AG8CuwLiIeBNA0s1Av2ZiKd72fHNxhYj4UNI9wMGSbgW+ApwF7EOWSD6aYugEPN7Cd78M6BIRs1qoZ2ZmZlWsmhO9hWkW7mOS9gUOAAZFxAJJ44DO6fbNwFHATOD2iAhlmdV1EfE/Rf0cCrTlBYQ7kc3EFcwvUe9msuTybaA+It5PMdwXEUe3YbyP0p+ZmZnVsKpfui3SA3gnJXlbAbvl7o0GDgWOZumM2wPAEZLWA5C0lqRNgCeBfSWtLakjcGSpASUdDhwE3NSK+MYBOwPfycXwBLCHpMIzhl0lNTd7aGZmZgbUXqJ3D7CKpKlkO1SfKNyIiHeAZ4FNImJiKnuW7Pm4sanNfcCGEfEaMJxsCfV+smfw8s4ovF4F+Cawf2GZtzkRsQQYQ/ZM4ZhU9iYwFLgpxfAEn3zez8zMzKxJPgKtBvkItKb5CDQzM6tQPgLNluq3xeY0NEwodxhmZma2gjnRq1CSbgc2LSr+UUTcW454zMzMbOXjRK9CRcTglmuZmZlZLau1zRhmZmZmNcObMWqQN2M0zZsxzMysQnkzhi3VuLiRxoHDyh3GSmdu/Yhyh2BmZtauvHRrZmZmVqWqMtGTtJGkOyU9L+kFSZdI6rSCx5yX/u0jaXqufE9JEyXNlDRL0qntMU6Je5tImpRe1vyMpJOXZywzMzOrbFWX6KWzYUcDd0REX6Af0A345XL22+ZlbkkbADcCJ0fEVsAewAmSVtSO2deA3dMZv7sCZ0vquYLGMjMzs5Vc1SV6wP7ABxHxJ/j4WLEzyBKsekkf70KQNE5Sf0mrSbo23X9K0iHp/lBJt0j6O9kxaN0kPSBpsqRphXrNOBUYFRGTUyxvAWcBP0z9j5J0RC6ewqxgW8ch9b84Ihalj6tSnf+/ZmZm1krVuBljW2BSviAi3pP0Mtn5sUcBP5e0IdAzIiZJOh94MCJOkLQGMFHS/an5IGCHiHg7zeoNTv2tAzwh6a4ovXV5W+C6orIGYJsWvsMHbRznY5J6A3cDWwA/jIhXW2pjZmZm1akaZ3wENJUQCRgHHJk+HwXckq4PIlvmnJLqdAY2Tvfui4i3c32cL2kqcD/QC1h/GWJpzXdoyzgfi4i5EbEDWaJ3vKRWtTMzM7PqU42J3jPAgHyBpO5Ab6Ae+LekHYAhwF8LVYDDI6Iu/W0cETPSvfm5ro4B1gX6p+fgXidLClsdC9CfbFYP4EPS/0F6trCwYaSt43xKmsl7BtirLe3MzMyselRjovcA0FXScQCSOgC/IXtWbgFZcncW0CMipqU29wLDUrKFpJ1K9N0DeCMiGiXtB2zSQiyXAUMl1aV+1ybbFHJeuj+HLPEDOATouIzjkPrfSFKXdL0m2eaPWa1pa2ZmZtWn6hK99BzbYOBISc8Dz5E983ZOqnIr8HXgb7lm55ElWVPTq1HOo2k3AAMkNZDNus1sIZbXgG8CV0qaBbwKXBoRD6cqVwH7SJpItku2MHvYpnFytgaelPQ08DBwcS6ZNTMzsxrjI9A+Q+kdeicDe0fEO+WKo0vXbtHr2GvKNfxKq2P9CGZMnlDuMMzMzNqq5BFoTvRqkM+6bZrPujUzswrls26riaTtgT8XFS+KiF1b077fFpvT0OCZKzMzs2rnRK8Cpefu6sodh5mZma3cqm4zhpmZmZll/IxeDfIzek3zM3pmZlah/IyeLdW4uJHGgcPKHcZKZ279iHKHYGZm1q68dGtmZmZWpZzomZmZmVWpFhM9SSHpN7nPZ0oa3kKbr0k6u4U6+0oaU+LeHEnrtBRbM30Pl3TmsrZfnn7T7zNT0nRJTxeOYmvHGFaVdL+kKZKGtGffZmZmVl1aM6O3CDisLYlXRNwVERcue1jLTlLZnjuUdDJwILBLRGwH7E0TD0im83eX1U5Ax4ioi4ibWxnX8oxnZmZmFao1id6HwJXAGcU3JK0r6TZJ9elvj1Q+VNLIdL25pCfS/XMlzct10U3SrWkG7AZJ+aToh5Impr8tUl+bSHpA0tT078apfJSk30p6CPhVar+NpHGSXpR0Wi7m76fZtumSTm9F+Y8lzZJ0P7BlC7/VOcB/RcR7ABHxbkRcl/qZI+lnkiaQncP7nfSbPJ1+w66SOqR4JWkNSR9J2ju1Hy9pF+AvQF2a0dtc0uclPSVpmqRrJa3a1HgtxG1mZmZVqLXP6F0GHCOpR1H5JcDvImIgcDhwdRNtLwEuSXVeLbq3E3A6sA2wGbBH7t57EbELMBL4fSobCVwfETsANwCX5ur3Aw6IiB+kz1sBXwB2AX4uqaOk/sC3gF2B3YDvSNqphfKvpzgPAwaW+oEkrQ6sHhEvlKoDfBARe0bEX4HRETEwInYEZgDfjoglwHPp99gTmATslZK3jSJiInAiMD4i6oBXgFHAkIjYnmwX9SklxjMzM7Ma06pEL81QXQ+cVnTrAGCkpCnAXUD3lPDkDQJuSdc3Ft2bGBH/ioiPgClAn9y9m3L/Dsr1Vejjz2TJUMEtKVEquDsiFkXEW8AbwPqp/u0RMT8i5gGjgb2aKd8rlS9Iv8FdTfw8BQJaeilhfql1uzRLNw04Bii82G482ZLv3sAFKbaBQH0T/W0JzI6I59Ln61K7psYzMzOzGtOWXbe/B74NrFbUflB6XqwuInpFxPtt6HNR7noJn3yvX5S4pkT5/Fb0XeqFgiVfNNjM2J+slCWC8yVt1ky1fIyjgO+lmbhfAJ1T+XiyBHMX4B/AGsC+wCNtjLt4PDMzM6sxrU70IuJt4G9kyV7BWOB7hQ+Smjp/9QmyZV3IlkFba0ju38fT9WO5Po4BJrShP8iSpUPT83CrAYPJEqvmygdL6pJmKg9uof8LgMskdQeQ1F3SSSXqrg68Jqlj+i4FTwK7Ax9FxAdkM53fTfEUmwn0KTzDCBwLPNxCjGZmZlYj2rpD9TfkEjuypdzLJE1NfT0CnFzU5nTgL5J+ANwNvNvKsVaV9CRZMnp0brxrJf0QeJPsubpWi4jJkkYBE1PR1RHxFGQbOkqU30yWbL1E08lW3h+AbkC9pEagkew3a8pPyZK6l4BpZIkfEbFI0lyyBJk05tGpTvH3+UDSt4Bb0m7jeuCKFmI0MzOzGrHCz7qV1BVYGBEh6evA0RFxyAod1Jrls26b5rNuzcysQpX1rNv+ZBs2BPwHOOEzGNOa0W+LzWloaOuqt5mZmVWaFZ7oRcR4YMcVPc5nSdJlfPJVMJC9QuZP5YjHzMzMrCllO0WikkXEqeWOwczMzKwlbXm9ipmZmZlVkBW+GcNWPt6M4Y0XZmZWVcq6GcNWMo2LG2kcOKzcYZTV3PoR5Q7BzMxshfPS7XKQNK/o81BJI1to83EdSetKelLSU5L2kjRH0jRJU9K/Lb6GRtI5ues+kqYv6/cxMzOz6uJEr7w+D8yMiJ3S7mSA/SKiDjgCuLQVfZzTchUzMzOrRU70VhBJB+dm6+6XtH7R/Trg18CX0wxel6IuugPv5OrfIWmSpGcKx6pJuhDoktrfkKp2kHRVqje2iX7NzMysRjjRWz6FJGuKpCnAubl7E4DdImIn4K/AWfmGETEF+Blwc0TURcTCdOuhtPz6MPCTXJMTIqI/MAA4TdLaEXE22akjdRFROC+3L3BZRGxL9oLqwzEzM7Oa5M0Yy2dhWmYFsufvyBIxgI2AmyVtCHQCZreyz/0i4i1JmwMPSBoXEfPIkrvBqU5vsoTu3020n52SSIBJQJ+2fCEzMzOrHp7RW3FGACMjYnvgu0DntjSOiBeA14FtJO0LHAAMiogdgaea6W9R7noJTubNzMxqlhO9FacH8Eq6Pr6tjSWtB2wKvJT6eiciFkjaCtgtV7VRUsflDdbMzMyqjxO9FWc4cIuk8cBbbWj3UHre7yHg7Ih4HbgHWEXSVOA84Ilc/SuBqbnNGGZmZmaAT8aoSV26dotex15T7jDKqmP9CGZMnlDuMMzMzNqDT8awpTp26kjHGj8ZonfP9codgpmZ2QrnRK8G9dticxoaPJtlZmZW7fyMnpmZmVmVcqJnZmZmVqWc6JmZmZlVKe+6rUHd11gzem22bbnD+Mz17rkeY8eMLncYZmZm7c27bm2pxsWNNA4cVu4wPnNza3ynsZmZ1R4v3ZqZmZlVqZpN9CRtIOmvkl6Q9Kykf0jqtwz9DJXUcxnaDZd0Zu7zKpLeknRBUb2rJW3Tyj7Pl/Sr3OdNJL0oaY22xmdmZmaVryYTPUkCbgfGRcTmEbENcA6w/jJ0NxRoMtGT1KEN/RwEzAKOSvEBEBEnRsSzrez7POAQSVunz5cAP42I/7QhDjMzM6sSNZnoAfsBjRFxRaEgIqZExHhJP5RUL2mqpF8ASOojaYakqyQ9I2mspC6SjgAGADdImpLK5kj6maQJwJGSvpP6e1rSbZK6lojpaLLE7GVgt0KhpHGSBqTreZLOlfQkMKi4g4hYCHwfuFzSl4DVI8Jn4JqZmdWoWk30tgMmFRdKOgjoC+wC1AH9Je2dbvcFLouIbfn/7d19tF1Vfe7x7yMEIQkQC2JJCIa3lALGAImIIOWlvsJozAUER0BS6b2lZVDlVhHp6LCVoa1VSwUpCvKiFgHFXMEMFTCAhKslOUkgCYSIrQGiUcTwjgkn4ekfax6zOZ73k312ztrPZ4w9stdca835WyvrnPMbc625JjwNnGz7ZqADmGN7ekm0ADbYPtr2jcA82zNtvxFYBZzdQ7s7AScA84EbqJK+nowDVto+wnaPU1vY/i6wHvgq8Nf9nYiIiIior3ZN9Hrz9vJZBiwFDqRK8AB+Zvv+8n0JMKWPem5q+H6IpIWSVgBzgJ7ea3IScJftF4FvAbN7uTW7uazvz+XAYturB7BtRERE1FS7vl7lQeCUHsoF/JPtL72iUJoCbGwo2gzs1Ef9LzR8vw54j+0HJM0Fju1h+/cBR0laU5Z3o7q9/INu222wvbmPdru8XD4RERHRxtq1R+9O4NWS/ndXgaSZwLPABySNL2WTJO3RT13PATv3sX5nYJ2kMVQ9eq8gaRfgaGBv21NsTwHOpffbtxERERED0pY9erYtaTbwb5IuBDYAa4APUT1/9+My8PV54AyqHrzeXAd8UdJv6WGABPD3wH3Ao8AKfj8p/F/AnbYbewxvAf5F0qsHd2QRERERW2QKtDaUKdAiIiJqJVOgxRZT99+Pjo4eB+1GREREjSTRG6Uk/T9gn27FH7V9WyviiYiIiG1PEr1RyvbsVscQERER27Z2HXUbERERUXsZjNGGMhgjIiKiVjIYI7bofKmTzpnntTqMEff44staHUJERMSIyq3bIZI0WdJdklZJelDSBwe5/92SZpTvayStkHR/+bxF0hRJK3vZ91WSLpW0suy3WNI+vdU1/KONiIiI0Sg9ekO3Cfhb20sl7QwskXSH7YeGWN9xtp/sWijTrv0eSdsDpwITgWm2X5a0F6+cdu0VdUVERER7SqI3RLbXAevK9+ckrQImSfp3qpkwjgMmAGfbXihpJ+Ba4CBgFX3PlfsKZY7cE4EdgXHAfGCd7ZdL+2u31nFFREREfSTR2wpK79uhVAkewPa23yTp3cDHgT8F/gp40fY0SdOApd2quUvSZmCj7SN6aOZIqh689aUH715JbwUWAP9he9kg6oqIiIg2kERvmCSNB74FfMj2s2WO3K6hnUuAKeX7McClALaXS1rerar+brfeYXt92X+tpD8Cji+fBZJOtb1ggHVFREREG0iiNwySxlAledfbbnxvx8by72ZeeY6H8y6bxmfwsL0R+B7wPUm/At5D1bsXERERAWTU7ZCp6rq7Glhl+18HsMs9wJyy7yHAtGG0fZikieX7q0pdjw61voiIiKin9OgN3VHAmcAKSfeXsov62P4K4Npyy/Z+YNEw2t4DuErSq8vyIuALw6gvIiIiaigzY7ShncaO96Qzr251GCNuzOLLWLX03laHERERsbVlZozYYswOYxjThrNETJ64R6tDiIiIGFFJ9NrQ1P33o6MjPVsRERF1l8EYERERETWVRC8iIiKippLoRURERNRURt22oV0mvMaT9j241WGMqMkT9+D2+fP63zAiImL0yajb2KLzpU46Z57X6jBG1ONtOMo4IiIit263EknXSHpC0sp+tjtW0lsalv9B0s8l3V8+/1zK75Y0o5c6TpK0TNIDkh6S9Jd91RURERHtKT16W891VLNTfLWf7Y4Fngd+1FB2ie3PDqSRMhvGlcCbbK8ty1OGUldERETUW3r0thLb9wDrG8sk/U3pcVsu6UZJU4BzgPNLj9tbB1K3pOclfULSfcARVAn6b0q7G22v3prHEhEREfWQRK+5LgQOtT0NOMf2GuCLVL1u020vLNud33C79R091DMOWGn7iJJQ3go8KukGSXMkNf4/9ldXREREtIkkes21HLhe0hnApj6260r8ptu+rYf1m4FvdS3Y/gvgBGAR8GHgmkHUFREREW0iiV5znQhcDhwOLJE01GciN9je3Fhge4XtS4C3AScPL8yIiIiooyR6TVJup062fRdwATABGA88B+w8jHrHSzq2oWg68OgwQo2IiIiayqjbrUTSDVQjaneXtBa4GDhT0q5ULzK8xPbTkr4D3CxpFjCUl9kJuEDSl4DfAi8Ac7fCIURERETNZGaMNrTT2PGedObVrQ5jRI1ZfBmrlt7b6jAiIiKaITNjxBZjdhjDmDabKWLyxD1aHUJERMSIS6LXhqbuvx8dHendioiIqLsMxoiIiIioqSR6ERERETWVRC8iIiKipjLqtg3tMuE1nrTvwa0Oo6kmT9yD2+fPa3UYERERIyGjbmOLzpc66Zw5lFf4jR6Pt9mo4oiIiJ7k1m1ERERETdUu0ZNkSV9rWN5e0q8lzS/Lr5M0X9IDkh6S9N1Sfq6k+xs+K0tdfzzEOL4racLWOSqQdKykZyQtk/SwpM82rJtbYj2hoWx2KTtla8UQERERo0vtEj2qKcEOkbRTWX4b8POG9Z8A7rD9RtsHARcC2L7c9vSuD3ArcL3tVUMJwva7bT899MPo0ULbhwKHAidJOqph3QrgfQ3LpwMPbOX2IyIiYhSpY6IH8D3gxPL9fcANDev2BNZ2Ldhe3n1nSccA7wX+uizvKOlaSStKj9pxpXyupHmSvi/pEUn/0lDHGkm7S5oiaZWkqyQ9KOn2riRU0kxJyyX9WNJnJK0cyMHZ/i1wPzCpoXgh8CZJYySNB/Yv20RERESbqmuidyNwuqQdgWnAfQ3rLgeulnSXpL+TNLFxx3K79VrgLNvPluJzAWy/gSpx/EqpG2A6cBrwBuA0SZN7iOcA4HLbBwNPAyeX8muBc2wfCWwe6MFJek2p856GYgM/AN4BzKLqkYyIiIg2VstEr/TSTaFKyr7bbd1twL7AVcCBwDJJr23Y5ArgP2z//4ayo4Gvlf0fBh4FppZ1C2w/Y3sD8BDw+h5C+pntrt61JcCUklDubPtHpfzrAzi0t0paDvwSmG/7l93W30h1y/Z0XtmLGREREW2ololecSvwWXpIeGyvt/1122cCi4FjACSdRZUgXtxtl17fTwNsbPi+mZ5fWdPTNn3V2ZuFtqdR9R7+laTpjSttLwIOAXa3/ZMh1B8RERE1UudE7xrgE7ZXNBZKOl7S2PJ9Z2A/4DFJ+wKfBObY3tStrnuAOWWfqcDewOrhBGf7KeA5SW8uRacPYt+fAP8EfLSH1R8DLhpObBEREVEPtX1hsu21wOd7WHU48AVJm6gS3S/bXizpS8A4YJ70is62bckTxQAAC3JJREFU84B/B74oaQWwCZhre2O37YbibOAqSS8AdwPPDGLfLwIflrRPY6Ht7w03qIiIiKiHTIHWQpLG236+fL8Q2NP2B5vdbqZAi4iIqJVMgbaNOlHSx6j+Hx4F5o5Eo1P334+OjntHoqmIiIhooSR6LWT7JuCmxjJJ7wA+3W3Tn9mePWKBRURERC0k0dvGlNe/3NbqOCIiImL0q/Oo24iIiIi2lsEYbSiDMSIiImolgzFii86XOumceV6rw2iqxxdf1uoQIiIiWi63biMiIiJqqpaJnqS9JN0i6RFJ/yXp85J2aHKbXe/DmyJpZUP50ZIWSXpY0mpJ526NdvpY/2lJK8vntOG0FREREaNb7RI9VdNVzAO+bfsAYCownmp6s+HUO+jb3JL+EPg6cI7tA4GjgA9IasqrUiSdCBwGTAeOAD4iaZdmtBURERHbvtolesDxwAbb1wLY3gycT5VgLZb0u1EIku6WdLikcZKuKeuXSZpV1s+V9E1J3wFulzRe0gJJSyWt6NquD+cC19leWmJ5ErgA+Eip/zpJpzTE09UrONh2uhwE/ND2JtsvAA8A7xzgvhEREVEzdUz0DgaWNBbYfhZ4DJgPvBdA0p7ARNtLgL8D7rQ9EzgO+IykcWX3I4GzbB8PbABm2z6sbPc59T3h7e/FAnRQJWR9GWw7XR4A3iVprKTdy76TB7BfRERE1FAdR90K6OmdMQLuBq4APk6V8H2zrHs78GeSPlyWdwT2Lt/vsL2+oY5PSToGeBmYBLwO+OUgYxnIMQymHQBs3y5pJvAj4NfAj4FNQ2g/IiIiaqCOPXoPAjMaC8pzapOBxcBvJE0DTgNu7NoEONn29PLZ2/aqsu6FhqrmAK8FDrc9HfgVVVI44FiAw6l69aBKwl5VYhTQNWBksO38ju1PlmN4WzmuRwayX0RERNRPHRO9BcBYSe8HkLQd8DmqZ+VepEruLgB2tb2i7HMbcF7X7VFJh/ZS967AE7Y7JR0HvL6fWC4H5kqaXurdjWpQyMVl/RqqxA9gFjBmiO3QdaylDUoyOw24fSD7RkRERP3ULtFzNdXHbOBUSY8AP6F65u2issnNwOnANxp2u5gqyVpeXo1yMT27HpghqYOq1+3hfmJZB5wBXClpNfAL4FLbPyybXAX8iaRFVKNku3oPB9VOgzHAQkkPAVcCZ9jOrduIiIg2lSnQRlB5h945wDG2n2pVHDuNHe9JZ17dquZHxJjFl7Fq6b2tDiMiImIk9DpgM4leG8pctxEREbWSuW7rRNIbgK91K95o+4iB7D91//3o6EhvV0RERN0l0RuFyiCS6a2OIyIiIrZttRuMERERERGVPKPXhvKMXkRERK3kGb3YovOlTjpnntfqMJrq8cWXtTqEiIiIlsut24iIiIiaSqIXERERUVOjOtGTtJekWyQ9Ium/JH1e0g797zmsNp8v/04ps2h0lR8taZGkhyWtLi9HHnY7faz/vqSnJc3vVr6PpPvKObmp2ecjIiIitl2jNtEr89LOA75t+wBgKjCeai7Z4dQ76OcWJf0h8HXgHNsHAkcBH5A0ezix9OMzwJk9lH8auKSck6eAs5sYQ0RERGzDRm2iBxwPbLB9LYDtzcD5VAnWYkm/G1Yq6W5Jh0saJ+masn6ZpFll/VxJ35T0HeB2SeMlLZC0VNKKru36cC5wne2lJZYngQuAj5T6r5N0SkM8Xb2Cg23nd2wvAJ5rLCvJ7/FU8/kCfAV4z0DrjIiIiHoZzaNuDwaWNBbYflbSY8B84L3AxyXtCUy0vUTSp4A7bX9A0gRgkaQflN2PBKbZXl969WaX+nYH/lPSre79XTQHUyVVjTqAg/o5hg2DbKc/uwFP295UltcCk4ZYV0RERIxyo7lHT0BPCZGAu4FTy/J7gW+W728HLpR0f9lmR2Dvsu4O2+sb6viUpOXAD6iSpdcNIZaBHMNg2hlIfd3lRYkRERFtajT36D0InNxYIGkXYDKwGPiNpGnAacBfdm0CnGx7dbf9jgBeaCiaA7wWONx2p6Q1VElhX7HMAG5tKDucqlcPYBMlqS63V7sGSAy2nf48CUyQtH3p1dsL+MUw6ouIiIhRbDT36C0Axkp6P4Ck7YDPUT0r9yJwI9VzcruWuWEBbgPOK8kWkg7tpe5dgSdK8nUc8Pp+YrkcmCtpeql3N6pBIReX9WuoEj+AWcCYIbbTp3LL9y6g63nAs4BbhlNnREREjF6jNtErSc1s4FRJjwA/oXrm7aKyyc3A6cA3Gna7mCrJWl5ejXIxPbsemCGpg6rX7eF+YlkHnAFcKWk1VS/apbZ/WDa5CvgTSYuAxt7DQbXTSNJCqlvSJ0haK+kdZdVHgf8r6adUz+xdPdA6IyIiol4y120TlHfonQMcY/upVsfTXea6jYiIqJVe57pNoteGZsyY4Y6Ojv43jIiIiNGg10RvNA/GqD1JbwC+1q14o+0jWhFPREREjC5J9LZhZRDJ9FbHEREREaNTbt22IUnPAav73TCaZXeqV+HEyMu5b62c/9bK+W+dZp/7J22/s6cV6dFrT6ttz2h1EO1KUkfOf2vk3LdWzn9r5fy3TivP/ah9vUpERERE9C2JXkRERERNJdFrT1e2OoA2l/PfOjn3rZXz31o5/63TsnOfwRgRERERNZUevYiIiIiaSqLXZiS9U9JqST+VdGGr46kzSZMl3SVplaQHJX2wlP+BpDskPVL+fU2rY60zSdtJWiZpflneR9J95fzfJGmHVsdYR5ImSLpZ0sPlZ+DIXPsjR9L55ffOSkk3SNox137zSLpG0hOSVjaU9Xi9q3Jp+Tu8XNJhzYwtiV4bkbQdcDnwLuAg4H2SDmptVLW2Cfhb238MvBk4t5zvC4EFtg8AFpTlaJ4PAqsalj8NXFLO/1PA2S2Jqv4+D3zf9oHAG6n+D3LtjwBJk4C/AWbYPgTYDjidXPvNdB3Q/T12vV3v7wIOKJ//A1zRzMCS6LWXNwE/tf3ftl8CbgRmtTim2rK9zvbS8v05qj90k6jO+VfKZl8B3tOaCOtP0l7AicCXy7KA44GbyyY5/00gaRfgGOBqANsv2X6aXPsjaXtgJ0nbA2OBdeTabxrb9wDruxX3dr3PAr7qyn8CEyTt2azYkui1l0nA4w3La0tZNJmkKcChwH3A62yvgyoZBPZoXWS192/ABcDLZXk34Gnbm8pyfgaaY1/g18C15bb5lyWNI9f+iLD9c+CzwGNUCd4zwBJy7Y+03q73Ef1bnESvvaiHsgy7bjJJ44FvAR+y/Wyr42kXkk4CnrC9pLG4h03zM7D1bQ8cBlxh+1DgBXKbdsSUZ8FmAfsAE4FxVLcLu8u13xoj+nsoiV57WQtMbljeC/hFi2JpC5LGUCV519ueV4p/1dVNX/59olXx1dxRwJ9JWkP1mMLxVD18E8rtLMjPQLOsBdbavq8s30yV+OXaHxl/CvzM9q9tdwLzgLeQa3+k9Xa9j+jf4iR67WUxcEAZebUD1cO5t7Y4ptoqz4NdDayy/a8Nq24FzirfzwJuGenY2oHtj9ney/YUqmv9TttzgLuAU8pmOf9NYPuXwOOS/qgUnQA8RK79kfIY8GZJY8vvoa7zn2t/ZPV2vd8KvL+Mvn0z8EzXLd5myAuT24ykd1P1amwHXGP7ky0OqbYkHQ0sBFaw5Rmxi6ie0/sGsDfVL+RTbXd/iDe2IknHAh+2fZKkfal6+P4AWAacYXtjK+OrI0nTqQbB7AD8N/DnVJ0LufZHgKR/BE6jGv2/DPgLqufAcu03gaQbgGOB3YFfAR8Hvk0P13tJvr9ANUr3ReDPbXc0LbYkehERERH1lFu3ERERETWVRC8iIiKippLoRURERNRUEr2IiIiImkqiFxEREVFTSfQiIiIiaiqJXkRERERNJdGLiIiIqKn/AdnRZvXeD9WrAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 648x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from statlearning import plot_feature_importance\n",
    "\n",
    "plot_feature_importance(xbst, labels=predictors)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The XGBoost functionately extends well beyond the <TT>Scikit-Learn API</TT>. Below, we construct a pure XGBoost implementation to select the number of boosting iterations by cross-validation and early stopping. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Selected number of boosting iterations: 646\n",
      "RMSE (CV): 0.0569\n",
      "CPU times: user 2min 40s, sys: 1.54 s, total: 2min 41s\n",
      "Wall time: 21.1 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "dtrain = xgb.DMatrix(X_train, y_train) # we need to convert the data to the format required by xgboost\n",
    "dtest  = xgb.DMatrix(X_test)\n",
    "\n",
    "param = {\n",
    "    'max_depth': 2, \n",
    "    'learning_rate': 0.1, \n",
    "    'subsample': 0.8,\n",
    "    'objective': 'reg:squarederror',  \n",
    "     }\n",
    "\n",
    "cv = xgb.cv(param, dtrain, num_boost_round = 1500, nfold=10, early_stopping_rounds=50, verbose_eval=False)\n",
    "\n",
    "print(f'Selected number of boosting iterations: {cv.shape[0]}')\n",
    "print(f'RMSE (CV): {cv.iloc[-1,0]:.4f}');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAe4AAAE9CAYAAADNvYHXAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deZhldX3n8fe39qpe6aaAhgYBQVxQAdsFRIMYI8EtauKSzTFOemYSFU3GRJOM0cyTiUbjlnE0jBtRgjMRF0J8XIZFJBqgWQVBAUVoaKCbrZveq+o7f9xzu6urq27d7q57zz1136/nKe89S53z/d0u/NzfOb9zTmQmkiSpGnrKLkCSJDXP4JYkqUIMbkmSKsTgliSpQgxuSZIqxOCWJKlC+souoBlnnXVWfutb3yq7DEmS2iVmWlCJHveGDRvKLkGSpI5QieCWJEk1BrckSRVicEuSVCEGtyRJFWJwS5JUIQa3JEkVYnBLklQhBrckSRVicEuSVCFdF9xX3LKOf/63n5VdhiRJ+6XrgvuLl93B+y+4ruwyJEnaL10X3IN9PWzfOVF2GZIk7ZeuC+6B/l62j42XXYYkSful64J7sL+XHTsNbklSNXVdcA94qFySVGEtC+6I+FxEPBgRN0+atywivhsRtxevB7Vq/zMZ6Oth5/gEmdnuXUuSdMBa2eP+AnDWlHnvBi7JzOOBS4rpthrs7wVgx5i9bklS9bQsuDPzCuDhKbNfBZxXvD8P+LVW7X8m9eDe7nluSVIFtfsc96GZuQ6geD2kzftnoK/W5O32uCVJFdSxg9MiYnVErImINevXr5+z7e46VG6PW5JUQe0O7gciYgVA8frgTCtm5rmZuSozV42Ojs5ZAbt63Aa3JKmC2h3cFwFvKt6/CfhGm/fPQP0ct4fKJUkV1MrLwS4AfgicEBFrI+ItwAeAl0TE7cBLium2Gix63Du9e5okqYL6WrXhzHzjDIte3Kp9NmP3qHJ73JKk6unYwWmtMuDlYJKkCuu64K4fKvcGLJKkKuq64B7os8ctSaqurgvuwX4vB5MkVVfXBfeA9yqXJFVY1wX37nPc9rglSdXTfcHt5WCSpArruuD2cjBJUpV1X3B7OZgkqcK6LrgHvRxMklRhXRfcfb1BhMEtSaqmrgvuiGCwr9dD5ZKkSuq64IbaTVi8HEySVEVdGdwD/b1eDiZJqqTuDO6+HnZ4jluSVEFdGdyDfb1s9xy3JKmCujK4a4fK7XFLkqqnK4N7sN9D5ZKkaurO4PZQuSSporoyuAe8HEySVFHdGdx9vezwcjBJUgV1ZXAP9vWw3R63JKmCujK4vQGLJKmqujK4B/t7HVUuSaqk7gxuD5VLkiqqK4PbQ+WSpKrqzuDu62GnPW5JUgV1ZXAP2uOWJFVUVwb3QHGOOzPLLkWSpH3SlcE92N9LJoyNG9ySpGrp2uAGHFkuSaqcrgzugb5as320pySparo6uL1fuSSparoyuIcGPFQuSaqmrgzu+jnubTsMbklStXRlcA/VB6d5jluSVDFdGdzDA30AbLXHLUmqmK4M7sHiHPc2e9ySpIrpyuAergf3jrGSK5Ekad+UEtwR8c6IuCUibo6ICyJiqJ37H3JwmiSpotoe3BFxBPB2YFVmngj0Am9oZw1DxTluB6dJkqqmrEPlfcBwRPQBI8B97dx5vcft4DRJUtW0Pbgz817gw8DdwDrgscz8TjtrcHCaJKmqyjhUfhDwKuAY4HBgQUT89jTrrY6INRGxZv369XNaw+7BaQa3JKlayjhU/svAzzNzfWbuBL4KnDZ1pcw8NzNXZeaq0dHROS3AwWmSpKoqI7jvBp4XESMREcCLgVvbWUBfbw99veHlYJKkyinjHPdVwFeA64AfFTWc2+46hvp7PcctSaqcvjJ2mpl/CfxlGfuuGxroc1S5JKlyuvLOaVDrcXsdtySparo2uIcHeh2cJkmqnK4N7sGBXgenSZIqp2uD28FpkqQq6t7gHuh1cJokqXK6N7gdnCZJqqDuDW4vB5MkVdB+BXfxVK9Ks8ctSaqiGYM7Iq6c9P6LUxZf3bKK2sTLwSRJVdSox71g0vunTVkWLailrQYHetnq5WCSpIppFNy5n8sqYXigz0PlkqTKaXSuemlEvJpauC+NiNcU8wNY0vLKWmyw38vBJEnV0yi4vwe8ctL7V0xadkXLKmqTof5exieSsfEJ+nq7dnC9JKliZgzuzHxzOwtpt+GBXgC27hhn0bDBLUmqhkajyl8REU+YNP3eiLgxIi6KiGPaU17rDBbB7W1PJUlV0qir+dfAeoCIeDnw28DvARcBn259aa1V73H7oBFJUpU0HFWemVuK968BPpuZ12bmZ4DR1pfWWkP99eC2xy1Jqo5GwR0RsTAieoAXA5dMWjbU2rJab2igdnrfS8IkSVXSaFT5x4AbgI3ArZm5BiAiTgbWtaG2lqr3uL0kTJJUJY1GlX8uIr4NHALcOGnR/UDlR5w7OE2SVEUzBndEnDJp8qSIve5yendLKmoTB6dJkqqo0aHyNcAtFCPL2fP+5Amc2aqi2sHBaZKkKmoU3H8MvBbYCnwZ+FpmPt6WqtpgaMBz3JKk6plxVHlmfjQzTwfeChwJXBIR/zciTmpbdS1kj1uSVEWz3uszM38OfAP4DvAc4EmtLqodhgdrBxt8tKckqUoaDU47FngD8CrgHmqHy/86M7e1qbaWWlAE95btBrckqToaneO+A7iJWm97I3AU8Af10eWZ+ZGWV9dCwwMGtySpehoF919RGz0OsLANtbRVT08wPNBrcEuSKqXRDVjeN9OyiFjQkmrabHiwjy3bHZwmSaqOhoPTIuKIiFgVEQPF9CER8T+A29tSXYstGOxjy/adZZchSVLTGj2P+x3U7lX+98C/R8SbgFuBYeBZ7SmvtexxS5KqptE57tXACZn5cEQcRW2w2gsz89/bU1rrjQz0eTmYJKlSGh0q35aZDwNk5t3AT+dTaAOMDPWx2cFpkqQKadTjXhkRn5g0fcjk6cx8e+vKao+RgV4e3byj7DIkSWpao+B+15Tpa1tZSBlGBvu59+EtZZchSVLTGl0Odl47CynDyGAvWz1ULkmqkFnvVT6fDQ96jluSVC1dHdwLBvvscUuSKqWrg7t+HXdmzr6yJEkdoNHgNAAiYhT4feDoyetn5u/t704jYinwGeBEavdD/73M/OH+bm9/LRjsYyKT7TvHGRqY9aOQJKl0zaTVN4DvA/8PmKvbjH0c+FZm/npxO9WROdruPqk/k3vLDoNbklQNzaTVSGb+6VztMCIWAy8E/gNAZu4ASrmYetczubeNsWzhYBklSJK0T5o5x31xRJw9h/s8FlgPfD4iro+Iz5T1tLGRXT1uB6hJkqqhmeA+h1p4b4uITcXPxgPYZx9wCvCpzDwZ2Ay8e+pKEbE6ItZExJr169cfwO5mNjzQC+DIcklSZcwa3Jm5KDN7MnOoeL8oMxcfwD7XAmsz86pi+ivUgnzqfs/NzFWZuWp0dPQAdjezBUP9AF7LLUmqjKZGZEXEK6mdlwa4PDMv3t8dZub9EXFPRJyQmT8BXgz8eH+3dyDscUuSqqaZy8E+ADwbOL+YdU5EnJ6Zex3e3gdvA84vRpT/DHjzAWxrvy0YqjXfHrckqSqa6XGfDZyUmRMAEXEecD3TnJduVmbeAKza39+fK8PFJWD2uCVJVdHsndOWTnq/pBWFlKE+qtwetySpKprpcf8NcH1EXAYEtXPd72lpVW1SD2573JKkqpg1uDPzgoi4nNp57gD+NDPvb3Vh7bDrOu7tc3VDOEmSWmvGQ+UR8eTi9RRgBbXLuO4BDi/mVd5AXw+9PeGhcklSZTTqcf8RsBr4u2mWJXBmSypqo4hgxEd7SpIqZMbgzszVxdtfzcxtk5dFxFBLq2qjkcE+e9ySpMpoZlT5D5qcV0kjg31sMbglSRUxY487Ig4DjgCGI+JkagPTABZT0mM4W2HhUB+bt+0suwxJkprS6Bz3S6k9enMl8JFJ8zcBf9bCmtpq4VA/m7Ya3JKkamh0jvs84LyIeG1mXtjGmtpq4XA/D23aNvuKkiR1gGau474wIl4GPA0YmjT/r1pZWLssHOrnrgc3lV2GJElNmXVwWkR8Gng9tQeDBPAbwBNaXFfbLBru5/GtDk6TJFVDM6PKT8vM3wUeycz3A6cCR7a2rPZZONzH4w5OkyRVRDPBvbV43RIRhwM7gWNaV1J7LRzq5/FtO8nMskuRJGlWzQT3xRGxFPgQcB1wF/DlVhbVTouG+8nEa7klSZXQzOC0/168vTAiLgaGMvOx1pbVPguH+wF4fNtOFgz1l1yNJEmNNboBy2saLCMzv9qaktprURHWm7aOcejSWVaWJKlkjXrcryheDwFOAy4tpl8EXA7Mi+De1eP2JiySpApodAOWNwMUh8efmpnriukVwCfbU17r1YN7kyPLJUkV0MzgtKProV14AHhSi+ppu4VDte8u9rglSVUw6+A04PKI+DZwAbXncL8BuKylVbXRwl3nuA1uSVLna2ZU+VuLgWovKGadm5lfa21Z7bNo0qhySZI6XTM97voI8nkxGG2q+jnuzdu8jluS1PkaXQ52ZWaeHhGbqB0i37UIyMxc3PLq2qB+jttD5ZKkKmg0qvz04nVR+8ppv96eHkYG+xycJkmqhEY97mWNfjEzH577csqxcKjPy8EkSZXQ6Bz3tdQOkcc0yxI4tiUVlaD2aE+DW5LU+RodKp83TwCbzUKDW5JUEU2NKo+Ig4DjgaH6vMy8olVFtdvCoX4PlUuSKmHW4I6I/wicA6wEbgCeB/wQOLO1pbXPwqF+Hnh06+wrSpJUsmZueXoO8GzgF5n5IuBkYH1Lq2qzhcN9bLbHLUmqgGaCe1tmbgOIiMHMvA04obVltdei4QE2eo5bklQBzZzjXhsRS4GvA9+NiEeA+1pbVnstGeln45YdZZchSdKsmrlX+auLt++LiMuAJcC3WlpVmy1ZMMDWHeNs3znOYH9v2eVIkjSjWQ+VR8THI+I0gMz8XmZelJnzqnu6ZGQAgMfsdUuSOlwz57ivA/4iIu6IiA9FxKpWF9VuSxcUwb3Z4JYkdbZZgzszz8vMs4HnAD8FPhgRt7e8sjZaUgT3xi0OUJMkdbZmetx1xwFPBo4GbmtJNSXxULkkqSqaOcdd72H/FXAL8KzMfEXLK2ujeo/7UQ+VS5I6XDOXg/0cODUzN8zljiOiF1gD3JuZL5/Lbe8rz3FLkqqimXPcn66HdkS8bw73fQ5w6xxub78t9lC5JKki9uUcN8Ar52KnEbESeBnwmbnY3oFaMNhHb094qFyS1PH2Nbinezb3/vgY8CfAxBxt74BEBEsWDHioXJLU8fY1uJ91oDuMiJcDD2bmtbOstzoi1kTEmvXrW/9Mk6UjAx4qlyR1vGZGlf9tRCyOiH5q9yrfEBG/fQD7fD7wyoi4C/gycGZEfGnqSpl5bmauysxVo6OjB7C75tjjliRVQTM97l/JzI3Ay4G1wJOAd+3vDjPzPZm5MjOPBt4AXJqZB/JFYE4sWTDgg0YkSR2vmeDuL17PBi7IzIdbWE9plowM8KjBLUnqcM1cx/0vEXEbsBX4g4gYBbbNxc4z83Lg8rnY1oFa6qFySVIFNHMd97uBU4FVmbkT2Ay8qtWFtdviEYNbktT5mhmc9hvAWGaOR8RfAF8CDm95ZW22ZGSAzdvH2DnWEVeoSZI0rWbOcf+3zNwUEacDLwXOAz7V2rLar36/ci8JkyR1smaCe7x4fRnwqcz8BjDQupLKsdQHjUiSKqCZ4L43Iv4BeB3wzYgYbPL3KuWghbXgfuTx7SVXIknSzJoJ4NcB3wbOysxHgWUcwHXcnWr5oiEAHt5kcEuSOlczo8q3AHcCL42ItwKHZOZ3Wl5Zmy1bNAjAQ5vm5Eo3SZJaoplR5ecA5wOHFD9fioi3tbqwdlu+K7jtcUuSOlczN2B5C/DczNwMEBEfBH4I/H0rC2u3JSMD9PaEwS1J6mjNnOMOdo8sp3g/V4/37BgRwbJFgx4qlyR1tGZ63J8HroqIrxXTvwZ8tnUllWf5oiEHp0mSOtqswZ2ZH4mIy4HTqfW035yZ17e6sDLY45YkdbqGwR0RPcBNmXkicF17SirP8kWD3LFuY9llSJI0o4bnuDNzArgxIo5qUz2l8lC5JKnTNXOOewVwS0RcTe3JYABk5itbVlVJlheHyjOTiHk3/k6SNA80E9zvb3kVHWL5okHGxpONW3eyZGTe3Y5dkjQPzBjcEXEccGhmfm/K/BcC97a6sDLU75728KbtBrckqSM1Osf9MWDTNPO3FMvmnfr9yh1ZLknqVI2C++jMvGnqzMxcAxzdsopKtHxSj1uSpE7UKLiHGiwbnutCOsGyXT1ug1uS1JkaBfc1EfH7U2dGxFuAa1tXUnkOXlzrca9/bGvJlUiSNL1Go8rfAXwtIn6L3UG9ChgAXt3qwsqwZGSAwf4eHnjU4JYkdaYZgzszHwBOi4gXAScWs/81My9tS2UliAgOXTpscEuSOlYz9yq/DLisDbV0hMOWjhjckqSO1cxjPbvKIfa4JUkdzOCe4tClw9xvcEuSOpTBPcVhBw3z8Kbt7BgbL7sUSZL2YnBPcejS2iXq6x/z7mmSpM5jcE9RD27Pc0uSOpHBPUU9uO9/xOCWJHUeg3uKw5aOAPa4JUmdyeCeYnRJ7X7lBrckqRMZ3FMM9veybOGgwS1J6kgG9zQOPWiYdY9sKbsMSZL2YnBPY+XyBdz70Oayy5AkaS8G9zRWHryAezYY3JKkzmNwT+PIgxewYeM2tu4YK7sUSZL2YHBP48iDFwJw70Oe55YkdRaDexpHLK9dy33PhsdLrkSSpD21Pbgj4siIuCwibo2IWyLinHbXMJt6j3ut57klSR2mr4R9jgF/nJnXRcQi4NqI+G5m/riEWqZ1+LIRInCAmiSp47S9x52Z6zLzuuL9JuBW4Ih219HIYH8vhy4dtsctSeo4pZ7jjoijgZOBq8qsYzpHHrzAc9ySpI5TWnBHxELgQuAdmblxmuWrI2JNRKxZv3592+tbuXyhPW5JUscpJbgjop9aaJ+fmV+dbp3MPDczV2XmqtHR0fYWCBw1uoC7NzzOxES2fd+SJM2kjFHlAXwWuDUzP9Lu/Tfr2MMWs33nBPc+bK9bktQ5yuhxPx/4HeDMiLih+Dm7hDoaeuKKxQDcsW6vo/iSJJWm7ZeDZeaVQLR7v/vquMNqwX3nuo286OmHl1yNJEk13jltBocvG2Gov5c7799UdimSJO1icM+gpyc49rBF3OmhcklSBzG4G3jiisXccb/BLUnqHAZ3A8cdtpi7HtjE+MRE2aVIkgQY3A09ccVidoxNeM9ySVLHMLgbOOGIJQDctvbRkiuRJKnG4G7gqUcdBMDNv3ik5EokSaoxuBtYMjLAUaMLuOVug1uS1BkM7lmceNQybja4JUkdwuCexdOecBA/vfcxtu8cL7sUSZIM7tmceNRBjE8kP7nXAWqSpPIZ3LN4mgPUJEkdxOCexXErFrNgsI/r7nyo7FIkSTK4Z9PX28Mpxx3M1bc/WHYpkiQZ3M147vGj3HTXw2zdMVZ2KZKkLmdwN+HZTxplbDy5/mceLpcklcvgbsKzjx8F4Jqfri+5EklStzO4m3DIkmGOOXQRP/yJ57klSeUyuJv0wqcdxpU/vp+xcR/xKUkqj8HdpDOfcTiPbt7BdXduKLsUSVIXM7ibdMbTVxABl9x0X9mlSJK6mMHdpOWLhjjl2IO59EaDW5JUHoN7H7z4mYdz9e3rWb9xW9mlSJK6lMG9D179vKMZn0guuuoXZZciSepSBvc+OPEJB3H84Yu58Ac/L7sUSVKXMrj3QUTw2tOO4cofP8D9j2wpuxxJUhcyuPfR604/lolMvnjZHWWXIknqQgb3PnrS4Ut40dNX8Jnv3ubNWCRJbWdw74f/dNZTuPehLfzLNXeXXYokqcsY3PvhV5+1kicetogPfOVGJiay7HIkSV3E4N4PvT09/MXrT+aWux/hn//tZ2WXI0nqIgb3fnrtqcfwjKOX8d/Ov5bHtuwouxxJUpcwuPdTT0/widWncv8jW/nzf7ym7HIkSV3C4D4Aq44b5ZxXPI0vXHo751/u5WGSpNYzuA/Qe99wCr904gredu4P+O4Na8suR5I0zxncB6i/r4cv/tEZPHnlUl73wUv5p+/Z85YktY7BPQeWLRzkm3/5Up57wiirP3klb/7497j3oc1llyVJmocis/OvQ161alWuWbOm7DJmNTY+wd9+9Sb+7us3kVl7mtjvnHk8pz35EAb6essuT5JUHTHjAoN77t314Cb+/uJbuOB7d7Jx604WDPbx7CeN8pSVSznhiCUcvmwBo0uGOGTJEAcvHmJ4oI+enhn/jSRJ3aezgjsizgI+DvQCn8nMDzRav2rBXbd5204u+9E6LrnxPq67cwO3rX2UzdvHpl13ZLCP4YFeFgz1MTLYz8hgL4N9vfT19TDQ10N/bw/99dfe2rw9lvX20NvbQ09Ab0/Q29NTvO7+6Sl+ps6vr9vTE/TEdMuDiCCi9oS0oHY53MyvQU8PBLXf6Zn8u5Oma69BT7GsZ8Z16vP3XnfX8llqiPCLkaRK6Zzgjohe4KfAS4C1wDXAGzPzxzP9TlWDe6qJiWTdI1u4/5GtrN+4lQcf3cZDm7axefsYW7eP7fG6ZfsYO8Ym2Fn/GZ+oTY/vnrdjbJyx8WTH2Dg7xiYYn0gmMqnAQZRS7P4SMPuXh576twFqLzHpCwJ7TO/eZm3d3dNRrD35C8/kLxFR/E99GVP2s2u9mFrD5O3u3ufudff8ksWU/exeN2Zp29QaZq6NmFTDjJ/B3l+kJn8mu/c/5TOasr/J8yfPnLps79/Zc7rRssa/M7W25ubvUXeTNe9Le6bOb6Y9U2tupj3M8G/TTHsa/05z/wZ7fjaN2zPdv8Hs7Zz5b61Rew47aJhVx43utZ8DsHfRhb653EuTngPckZk/A4iILwOvAmYM7vmipyc4YvkCjli+oKX7yUzGJ6b+1IMdJor3U38mprxO/r0EMim+GOSu9+yaN9NrbZ+ZtW1MTOx+Zcp0kkxMQNJoe/X97t7uru03UcNE0ZD68qk11WuYKL791NYp1q3PK+aTuz+X+uee7K5l97q795XTbre+7u7p+vv6yvX9TK1h8no5Zb1Jv75HbfV59c+GorZd29hru3vXtvd2p69h8mdQ38+en8HUz2vP/dbn7TE96Yvp3st277PR/MkTe//OlG1O/pUZlu09f8/pZtoz3e+oOl7+7KP48rvObMu+ygjuI4B7Jk2vBZ5bQh3zVkTQ1xs4Hk6aH2b7YjB1/p7Lmv8yM/PvNP7C0mjZ7lUO/MvMHvs7oHY2+wVsj19qWNvikQHapYzgnq77v9d3zIhYDawGOOqoo1pdkyR1rL0Poc94FFVdoIzruNcCR06aXgncN3WlzDw3M1dl5qrR0Tk9byBJUmWVEdzXAMdHxDERMQC8AbiohDokSaqcth8qz8yxiHgr8G1ql4N9LjNvaXcdkiRVURnnuMnMbwLfLGPfkiRVmfcqlySpQgxuSZIqxOCWJKlCDG5JkirE4JYkqUIMbkmSKqQSz+OOiPXAL+ZwkwcDG+Zwe1Vi27uTbe9Otr26NmTmWdMtqERwz7WIWJOZq8quowy23bZ3G9tu2+cbD5VLklQhBrckSRXSrcF9btkFlMi2dyfb3p1s+zzUlee4JUmqqm7tcUuSVEldF9wRcVZE/CQi7oiId5ddz1yLiM9FxIMRcfOkecsi4rsRcXvxelAxPyLiE8VncVNEnFJe5QcmIo6MiMsi4taIuCUizinmd0PbhyLi6oi4sWj7+4v5x0TEVUXb/09EDBTzB4vpO4rlR5dZ/1yIiN6IuD4iLi6mu6LtEXFXRPwoIm6IiDXFvHn/Nw8QEUsj4isRcVvx3/2p3dL2rgruiOgFPgn8KvBU4I0R8dRyq5pzXwCmXvv3buCSzDweuKSYhtrncHzxsxr4VJtqbIUx4I8z8ynA84A/LP5tu6Ht24EzM/OZwEnAWRHxPOCDwEeLtj8CvKVY/y3AI5l5HPDRYr2qOwe4ddJ0N7X9RZl50qRLn7rhbx7g48C3MvPJwDOp/ft3R9szs2t+gFOBb0+afg/wnrLrakE7jwZunjT9E2BF8X4F8JPi/T8Ab5xuvar/AN8AXtJtbQdGgOuA51K7+URfMX/X3z7wbeDU4n1fsV6UXfsBtHkltf+TPhO4GIguavtdwMFT5s37v3lgMfDzqf923dD2zOyuHjdwBHDPpOm1xbz57tDMXAdQvB5SzJ+Xn0dx+PNk4Cq6pO3FoeIbgAeB7wJ3Ao9m5lixyuT27Wp7sfwxYHl7K55THwP+BJgoppfTPW1P4DsRcW1ErC7mdcPf/LHAeuDzxSmSz0TEArqj7V0X3DHNvG4eVj/vPo+IWAhcCLwjMzc2WnWaeZVte2aOZ+ZJ1HqfzwGeMt1qxeu8aXtEvBx4MDOvnTx7mlXnXdsLz8/MU6gdCv7DiHhhg3XnU9v7gFOAT2XmycBmdh8Wn858anvXBfda4MhJ0yuB+0qqpZ0eiIgVAMXrg8X8efV5REQ/tdA+PzO/WszuirbXZeajwOXUzvMvjYi+YtHk9u1qe7F8CfBweyudM88HXhkRdwFfpna4/GN0R9vJzPuK1weBr1H70tYNf/NrgbWZeVUx/RVqQd4Nbe+64L4GOL4YcToAvAG4qOSa2uEi4E3F+zdRO/9bn/+7xYjL5wGP1Q8zVU1EBPBZ4NbM/MikRd3Q9tGIWFq8HwZ+mdpAncuAXy9Wm9r2+mfy68ClWZz4q5rMfE9mrszMo6n993xpZv4WXdD2iFgQEYvq74FfAW6mC/7mM/N+4J6IOKGY9WLgx3RB24HuGpxW/Pd5NvBTaucA/7zselrQvguAdcBOat8y30LtHN4lwO3F67Ji3aA2yv5O4EfAqrLrP4B2n07t0NdNwA3Fz9ld0vZnANcXbb8ZeG8x/1jgauAO4J+BwWL+UDF9R7H82LLbMEefwxnAxd3S9qKNNxY/t9T//6wb/uaL9pwErCn+7r8OHKdJvy4AAAQ1SURBVNQtbffOaZIkVUi3HSqXJKnSDG5JkirE4JYkqUIMbkmSKsTgliSpQgxuqQ0iYrx4gtONEXFdRJw2x9v/synTP5ij7a6KiE8U78+Yy7oj4uiI+M3p9iVpZl4OJrVBRDyemQuL9y8F/iwzf6kV22+ViHgf8Hhmfngffqcvd98zfOqyM4D/mpkvn5sKpe5gj1tqv8XUHjVZf07whyLi5uK5yq+fZf6KiLii6L3fHBEviIgPAMPFvPOL9R4vXs+IiMsnPbf4/OIuc0TE2cW8K4tnFV88tdDi9y8uHtzyn4F3Fvt5QXHHtgsj4pri5/nF77wvIs6NiO8A/1j0rL9fHGmYfLThA8ALiu29s76vYhvLIuLrUXt28r9HxDMmbftzRZt+FhFvL+YviIh/LY5o3Fz/vKT5qG/2VSTNgeGoPb1riNrjBs8s5r+G2h2gngkcDFwTEVcAp80w/zepPaLyr6P2fPmRzPx+RLw1aw8Zmc7JwNOo3Zv534DnR8Qaao86fGFm/jwiLmhUfGbeFRGfZlKPOyL+idozr6+MiKOoPTKz/nCTZwGnZ+bWiBgBXpKZ2yLieGp391tF7aEQu3rcRQ+87v3A9Zn5axFxJvCPxecB8GTgRcAi4CcR8Slqz6C/LzNfVmxrSaP2SFVmcEvtsbUerBFxKrWe6InUbtV6QWaOU3tAwveAZzeYfw3wuag9UOXrmXlDE/u+OjPXFvu+gdrz2h8HfpaZPy/WuQBYPf2vz+iXgacWHXiAxVHcOxu4KDO3Fu/7gf8ZEScB48CTmtj26cBrATLz0ohYPimM/zUztwPbI+JB4FBqt7H8cER8kNptT7+/j22RKsND5VKbZeYPqfWiR5n+cYPMND8zrwBeCNwLfDEifreJXW6f9H6c2hf2mfa7L3qAUzPzpOLniMzcVCzbPGm9dwIPUDt6sAoYaGLbjR7DuFd7MvOn1Hr5PwL+JiLeuw/tkCrF4JbaLCKeDPQCDwFXAK+PiN6IGKUWylfPND8inkDt+dP/m9rT0E4pNruz6IU36zbg2OLcNUAz54Q3UTs8Xfcd4K2T2jXTofolwLrMnAB+h1rbp9veZFcAv1Vs9wxgQzZ4vnpEHA5sycwvAR9m9+cizTseKpfao36OG2q9yTdl5nhEfA04ldoTnhL4k8y8v8H8NwHvioid1A5313vc5wI3RcR1WXusZUPFuec/AL4VERuofVmYzb8AX4mIVwFvA94OfDIibqL2/yVXUBvANtX/Ai6MiN+g9rjNem/8JmAsIm4EvkDtCWd17wM+X2x7C7sf1TiTpwMfiogJak/G+y9NtEeqJC8Hk7pURCzMzMeLUeafBG7PzI+WXZekxjxULnWv3y+OAtxC7XD2P5Rcj6Qm2OOWJKlC7HFLklQhBrckSRVicEuSVCEGtyRJFWJwS5JUIQa3JEkV8v8BUCpdL8JimJoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(8,5))\n",
    "plt.plot(cv.iloc[:,0])\n",
    "ax.set_ylabel('Cross-validation RMSE')\n",
    "ax.set_xlabel('Boosting iterations')\n",
    "sns.despine()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Additive Boosting\n",
    "\n",
    "This is an advanced specification. Since gradient boosting is an additive model fit by forward stagewise additive modelling, nothing stops us from fitting a gradient boosting model to the residuals of a linear regression specification, therefore boosting the linear model with additive trees. Below we fit the linear regression model using the Lasso.\n",
    "\n",
    "The only disadvantage is that there are no immediately available functions to add this model to our stack. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Best parameters found by randomised search: {'subsample': 0.8, 'n_estimators': 750, 'max_depth': 4, 'learning_rate': 0.01} \n",
      "\n",
      "CPU times: user 22.1 s, sys: 157 ms, total: 22.3 s\n",
      "Wall time: 2min 27s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "y_fit = lasso.predict(X_train)\n",
    "resid = y_train - y_fit\n",
    "\n",
    "model = xgb.XGBRegressor()\n",
    "\n",
    "tuning_parameters = {\n",
    "    'learning_rate': [0.01, 0.05, 0.1],\n",
    "    'n_estimators' : [250, 500, 750, 1000, 1500],\n",
    "    'max_depth' : [2, 3, 4],\n",
    "    'subsample' : [0.6, 0.8, 1.0],\n",
    "}\n",
    "\n",
    "gb_search = RandomizedSearchCV(model, tuning_parameters, n_iter = 16, cv = 5, return_train_score=False, n_jobs=4,\n",
    "                              random_state = 20)\n",
    "gb_search.fit(X_train, resid)\n",
    "\n",
    "abst = gb_search.best_estimator_\n",
    "\n",
    "\n",
    "print('Best parameters found by randomised search:', gb_search.best_params_, '\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model Stacking\n",
    "\n",
    "Model stacking is a learning method that aims to improve predictive accuracy by combining predictions from multiple models. (Note that the <TT>mlxtend</TT> package needs to be installed, see the beginning of the tutorial.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 4min 17s, sys: 2.8 s, total: 4min 20s\n",
      "Wall time: 48.2 s\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "StackingCVRegressor(cv=10,\n",
       "                    meta_regressor=LinearRegression(copy_X=True,\n",
       "                                                    fit_intercept=True,\n",
       "                                                    n_jobs=None,\n",
       "                                                    normalize=False),\n",
       "                    n_jobs=None, pre_dispatch='2*n_jobs', random_state=None,\n",
       "                    refit=True,\n",
       "                    regressors=[LinearRegression(copy_X=True,\n",
       "                                                 fit_intercept=True,\n",
       "                                                 n_jobs=None, normalize=False),\n",
       "                                Pipeline(memory=None,\n",
       "                                         steps=[('scaler',\n",
       "                                                 StandardScaler(copy=True,\n",
       "                                                                with_mean=True,...\n",
       "                                             min_child_weight=1, missing=nan,\n",
       "                                             monotone_constraints='()',\n",
       "                                             n_estimators=1000, n_jobs=8,\n",
       "                                             num_parallel_tree=1,\n",
       "                                             objective='reg:squarederror',\n",
       "                                             random_state=0, reg_alpha=0,\n",
       "                                             reg_lambda=1, scale_pos_weight=1,\n",
       "                                             subsample=0.6, tree_method='exact',\n",
       "                                             validate_parameters=1,\n",
       "                                             verbosity=None)],\n",
       "                    shuffle=True, store_train_meta_features=False,\n",
       "                    use_features_in_secondary=False, verbose=0)"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "models = [ols, lasso, ridge, xbst]\n",
    "\n",
    "stack = StackingCVRegressor(models, meta_regressor = LinearRegression(), cv = 10)\n",
    "stack.fit(X_train.values, y_train.ravel())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below we plot the coefficients of the OLS metamodel."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEJCAYAAABlmAtYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAWRklEQVR4nO3de5hkdX3n8fcHEOXmLTjK6AiKqAECgwKLq6IRMMQMriJ5QMkKu9HVrIurWXQ1ui66uUi8EEWDYlTEaCQqKBkjIjfxBqIwwMAutxUdGGQGUQRBRP3mj3Mayqa6fz1DV1dPz/v1POfpOr/zO6e+p6q6PnV+53R1qgpJkqazybgLkCTNf4aFJKnJsJAkNRkWkqQmw0KS1GRYSJKaDAvNS0kqySHjrmNdJFme5KR5UMcxSW7uH8Mjh7X18yvXYZvXJzl6ZEVr3ot/Z6EHon9zPGLIogurap8Zrr9tVS2b1P4Y4CdVdfds1DnN/V8PfKCq3j0L21oO3FJVRz7QbT2AGnYFLgcOBr4N3AbsOKRtU+DBVfXjGW73UcDPq+rOWaz1GOCQqtp1trap0dls3AVoQTgL+I+T2n75QDZYVT96IOtvxJ7U//xC9Z8Ek9yvrXfHTDdaVWtnqT5toByG0my4u6p+NGm6dWJhklcluTrJL5KsTfKVJJv1nyyPAP6oHx6pJM/t17l3GCrJDv38YUm+luSuJJck2S3Jrkm+leTnSb6R5AkD97tjki8m+VG//OIkywaWnwdsD7xr4v4Hlv37/r7uTHJjkhOSPHRg+ZZJTkpyRz+88xczeaCS7JPknL6e25KcnWRxv+zBSf6u394vklyQ5FmT1t85yZeS3J5kTZJ/6o/CJj6pn9Z3/U2/T/drm+g7eRgqyRFJLk9yd1/DSQPLfmsYKsnDkpzY13B7/1jtObD8yP6x2S/Jyn5/z514fvrhsf8N7DLw3B/ZLxv6epnJ46vRMSw0Uv0byAeBtwNPAfYHzugXvxv4Z7ojk+366VvTbO7twLHAHsBPgU8DxwNvAfYGHgK8f6D/1sCXgQOA3YHPA6cmeWq//GDgBuAdA/dPkt8DzgRO79c7GFgKfGxg2+/ut/sSYL++pn0bj8XuwLnAtcAzgX36/Z94I/xb4FDgP/fbuxw4I8lEXdsB5wMr+/3dv9/H05Ns0tf0yn5bE/szrG1Yba8CPgx8HNgNeAFwxRR9A3wJeCywrK/1fOCciVp7Dwbe3O/PM4CHAx/ql50CvAe4aqCuUxqvF41TVTk5rfcEnAT8im5IY3A6tl9+MN0Y+TbTrL98SHvRjWcD7NDPv2pg+bK+7eCBtiOBOxr1XgC8dWD+euDoSX1OBj46qW1pf3+L6N6g7wYOH1i+NV2AnTTNfX8KuGCKZVvRDd29fKBtU+A64C/7+XcAZ09a7xF9XXv384d0v9a/1WdY2zHAyoH5G4B3TlP7vY8T8Lz+Od5iUp8VwBsHnosCnjKw/PB+HzcZVsNMXi9O45s8tNNsOB/4L5Paftr//CrwA+D7Sb5C94n91Kq6fT3u57KB2zf3Py+f1LZVki2r6s4kW9ENdSyj++T6ILqjj8HtDPN04ElJDh1oS/9zR+BOYHO6k8UAVNUdSQZrGWYP7hsSmmzHvr5vDmzz10m+Dew8UNe+SYada9gR+E7j/odKsojuKOHsGa7ydGBLYG13kHGvh/R1TLi7qq4amF9Nt48PB25luNl8vWgWGRaaDXdW1bXDFlTV7UmeRjdEcwDdsMRfJ9mrqlav4/3cM7jpadomhlffDRwIHA1cQ/cmfzLdG/10NgH+AThuyLIb6YZH1kdmsGzY5YmD+/Uluv2Z7OYhbbNR1zCb9Pf37CHLfjZw+1eTlk1+fu5nll8vmkWes9DIVdWvquqcqnoz3Xj4VnSf9qEblth0RHf9LODkqvp8VV1GN9Sy46Q+w+7/YmCXqrp2yHQX3TmHe+jOOQDQH8W0LgG9mG4IZ5hr+1ruPaGdZFO6sf4rB+sCfjCkrvX+5F1VN9OF4H4zXOVi4NHAb4bUsWYd7nroc994vWhMDAvNhgcnecyk6VEASZYl+e9J9kiyPfAyYBvg//brXg/smuQpSbZN8qBZrOtq4MVJntaftP5HuqGSQdcDz07y2CTb9m3HAnsn+VBf95P6/fgwdENOwEeBY5MckGQXupPfrdB7F7BHfxXR7v0+vyLJ46vq58AJwDuTvCDJ7/bzjwb+vl//g8DD6E4E/7skT0yyf7+9bdb/YQLgr4DXJXl9kicnWZrkf0zR9yy64bIvJvnDJE9I8owkb08y7GhjKtcD2/fPz7b91WCt14vGxLDQbNgfuGnSdEm/7KfAi+jeYP4f3RDKK6rq6/3yj9C9EXwXWEt3ldBs+XNgDfB1uquiLuhvD3obsITuRPJagP4oZF+6E+tfAy4F/obfHuo5mu7KptP6nyvpzt1MqapW0D1WT+1ruRA4jPuG0v4n3dVRH6c7WbwbcGBV3dSvv5ru8fkN3RVCV9AFyN39tN6q6gTgNXRXTq3st7/LFH2L7mqpc+iev6v6up9Cd15ipj4P/CvduZK1wEtpv140Jv4FtySpySMLSVKTYSFJajIsJElNhoUkqWlB/lHegQceWGec4dfJSNI6mvIPNBfkkcUtt9wy7hIkaUFZkGEhSZpdhoUkqcmwkCQ1GRaSpCbDQpLUZFhIkpoMC0lSk2EhSWoyLCRJTYaFJKnJsJAkNS3ILxLUxun5yw5m1eo14y5DGqslixdx5vJTZ327hoUWjFWr13DPXkeNuwxprFZddPxItuswlCSpybCQJDUZFpKkJsNCktRkWEiSmgwLSVKTYSFJajIsJElNhoUkqcmwkCQ1GRaSpCbDQpLUZFhIkprGHhZJHpfki0muSXJdkvcl2TzJc5MsH9J/WZJLklya5MokrxpH3ZK0MRlrWCQJcCrwharaCXgysDXwV1P0fxBwInBQVe0O7AGcNzfVStLGa9z/z+J5wC+q6uMAVfXrJK8Hvg+cO6T/NnQ1/7jvfzdw1RzVKkkbrXEPQ+0CfG+woap+BvwQeNLkzlV1K3A68IMk/5Tk8CTj3gdJWvDG/UYboNahnap6BbAf8B3gaOBjI6tOkgSMPyyuAPYcbEjyUGAJcN1UK1XV5VV1HHAA8JKRVihJGntYnA1smeTlAEk2Bd4DnATcOblzkq2TPHegaSnwg9GXKUkbt7GGRVUV8GLgj5NcA1wN/AL4i77LfklumJjorn56Y5KrkqwA3g4cOYbSJWmjMu6roaiqVcBBQxadB2wxpP3rIy1IknQ/4x6GkiRtAAwLSVKTYSFJajIsJElNhoUkqcmwkCQ1GRaSpCbDQpLUZFhIkpoMC0lSk2EhSWoyLCRJTYaFJKlp7N86K82WJYsXseqi48ddhjRWSxYvGsl2DQstGGcuP3XcJUgLlsNQkqQmw0KS1GRYSJKaDAtJUpNhIUlqMiwkSU2GhSSpybCQJDUZFpKkJsNCktRkWEiSmgwLSVKTYSFJavJbZ7VgPH/ZwaxavWbcZUhzYsniRXP6TcuGhRaMVavXcM9eR427DGlOzPX/bnEYSpLUZFhIkpoMC0lSk2EhSWoyLCRJTYaFJKnJsJAkNRkWkqQmw0KS1GRYSJKaDAtJUpNhIUlqMiwkSU3ThkWSJUm+n+SR/fwj+vntk+yUZHmS65J8L8m5Sfbt+x2ZZG2SFUmuSPK5JFvOVtFJliZ5wWxtT5I0vWnDoqpWAScA7+yb3gmcCNwMfAk4sap2rKqnA0cBTxxY/ZSqWlpVuwC/BA6dxbqXAoaFJM2RmQxDHQfsk+R1wLOA9wCHA9+uqtMnOlXVyqo6afLKSTYDtgJ+0s9vn+TsJJf1Px/faP/jJCuTXJrk/CSbA+8ADu2PXGYzhCRJQzTDoqruAd5AFxqvq6pfArsAFzdWPTTJCuBG4JHAv/TtHwBOrqrdgE8B72+0vw34g6raHXhhf/9v474jl1NmtquSpPU10xPcfwjcBOw6bGGS0/pP/4P/4++UqloKPAa4nC5wAJ4BfLq//Um6o5Xp2r8JnJTklcCmM6xXkjSLmmGRZClwALAP8Pok2wFXAE+b6FNVLwaOpDuC+C1VVXRHFftOcRc1XXtVvRp4K7AEWJHkd1o1S5JmV+tqqNCd4H5dVf0QeBfwbrojgGcmeeFA9+mudnoWcF1/+1vAYf3tw4FvTNeeZMequrCq3gbcQhcatwPbNPdOkjQrNmssfyXww6r6aj//93RHEHsDy4D3Jvk7uqujbgf+cmDdQ5M8iy6QbujXA3gt8LEkbwDWAv+p0f6uJDsBAc4GLgV+CLypPyfyN563kKTRmjYsqupEuktlJ+Z/DTx9oMvQy1f7q6JOmmLZ9cDz1qH94CGbuRXYa6q6JUmzy7/gliQ1GRaSpCbDQpLUZFhIkpoMC0lSk2EhSWoyLCRJTYaFJKnJsJAkNRkWkqQmw0KS1GRYSJKaWt86K20wlixexKqLjh93GdKcWLJ40Zzen2GhBePM5ae2O0laLw5DSZKaDAtJUpNhIUlqMiwkSU2GhSSpybCQJDUZFpKkJsNCktRkWEiSmgwLSVKTYSFJajIsJElNhoUkqcmwkCQ1+RXlWjCev+xgVq1eM+4ypHstWbxowXx1vmGhBWPV6jXcs9dR4y5DutdC+mdcDkNJkpoMC0lSk2EhSWoyLCRJTYaFJKnJsJAkNRkWkqQmw0KS1GRYSJKaDAtJUpNhIUlqMiwkSU2GhSSpaaRhkeSOUW5fkjQ3PLKQJDXNeVgkOSjJhUkuSXJWkkf37c9JsqKfLkmyTZLtkpzft61M8uy+70uTXN63HTvX+yBJG5txHFl8A9inqvYAPgO8sW8/GnhNVS0Fng3cBbwM+ErftjuwIsli4FjgecBSYK8kL5rjfZCkjco4/lPe44BTkmwHbA58v2//JvDeJJ8CTq2qG5JcBHwsyYOAL1TViiTPA86rqrUAff99gS/M+Z5I0kZiHEcWxwMfqKrfA14FPASgqt4JvALYArggyVOr6ny6ILgR+GSSlwMZQ82StFEbx5HFw+je/AGOmGhMsmNVXQ5cnuQZwFOT3AXcWFUfSbIV8DS6Iaj3JdkW+AnwUroAkiSNyKjDYsskNwzMvxc4BvhskhuBC4An9Mtel+T3gV8DVwJfBg4D3pDkHuAO4OVVdVOSNwPn0h1l/GtVfXHE+yFJG7WRhkVVTTXMdb8396o6aki/T/TT5L6fBj79wKqTJM2Uf2chSWoyLCRJTYaFJKnJsJAkNRkWkqQmw0KS1GRYSJKaDAtJUpNhIUlqMiwkSU2GhSSpybCQJDWN4yvKpZFYsngRqy7y2+o1fyxZvGjcJcwaw0ILxpnLTx13CdKC5TCUJKnJsJAkNRkWkqQmw0KS1GRYSJKaDAtJUpNhIUlqMiwkSU2GhSSpybCQJDUZFpKkJsNCktRkWEiSmgwLSVKTX1E+xPOXHcyq1WvGXYbW0ZLFi/yacmlEDIshVq1ewz17HTXuMrSO/MdH0ug4DCVJajIsJElNhoUkqcmwkCQ1GRaSpCbDQpLUZFhIkpoMC0lSk2EhSWoyLCRJTYaFJKnJsJAkNRkWkqSmkYVFkl8nWZFkZZJ/SfLwvn1xks9Nsc55SfYcVU2SpPUzyiOLu6pqaVXtCtwKvAagqlZX1SEjvF9J0iybq2GobwOPBUiyQ5KV/e0tknwmyWVJTgG2mFghyZ8mubo/2vhIkg/07Y9K8vkkF/XTM+doHyRpozXyf36UZFNgP+CjQxb/GXBnVe2WZDfg4n6dxcD/Ap4G3A6cA1zar/M+4Liq+kaSxwNfAX53tHshSRu3UYbFFklWADsA3wO+OqTPvsD7AarqsiSX9e17A1+rqlsBknwWeHK/bH9g5yQT23hokm2q6vaR7IUkafTnLIDtgc3pz1kMUUPaMqRtwibAM/rzIUur6rEGhSSN1sjPWVTVbcBrgaOTPGjS4vOBwwGS7Ars1rd/B3hOkkck2Qx4ycA6ZwL/bWImydJR1S5J6szJCe6quoTunMNhkxadAGzdDz+9kS4kqKobgb8GLgTOAq4EbuvXeS2wZ39S/Erg1aPfA0nauI3snEVVbT1p/qCB2V37tru4f4BM+HRVndgfWZxGd0RBVd0CHDr7FUuSpjKf/4L7mP4E+Urg+8AXxlyPJG20Rn7p7PqqqqPHXYMkqTOfjywkSfOEYSFJajIsJElNhoUkqcmwkCQ1GRaSpCbDQpLUZFhIkpoMC0lSk2EhSWqat1/3MU5LFi9i1UXHj7sMraMlixeNuwRpwTIshjhz+anjLkGS5hWHoSRJTYaFJKnJsJAkNRkWkqQmw0KS1GRYSJKaDAtJUpNhIUlqMiwkSU2GhSSpybCQJDWlqsZdw6xLshb4wbjrmMK2wC3jLmKGNqRaYcOqd0OqFTaseq11/d1SVQcOW7Agw2I+S/Ldqtpz3HXMxIZUK2xY9W5ItcKGVa+1jobDUJKkJsNCktRkWMy9E8ddwDrYkGqFDaveDalW2LDqtdYR8JyFJKnJIwtJUpNhIUlqMixGLMkjk3w1yTX9z0cM6bM0ybeTXJHksiSHztda+35nJPlpkuVjqPHAJFcluTbJm4Ysf3CSU/rlFybZYa5rnFRPq959k1yc5FdJDhlHjQO1tGr98yRX9q/Rs5NsP446B+pp1fvqJJcnWZHkG0l2HkedfS3T1jrQ75AklWT+XU5bVU4jnIC/Bd7U334TcOyQPk8GdupvLwZuAh4+H2vtl+0HHAQsn+P6NgWuA54IbA5cCuw8qc9/BT7U3z4MOGWMz/1M6t0B2A04GThkntf6+8CW/e0/2wAe24cO3H4hcMZ8rbXvtw1wPnABsOe4HtupJo8sRu8/AJ/ob38CeNHkDlV1dVVd099eDawBHjVnFd6nWStAVZ0N3D5XRQ3YG7i2qv5/Vf0S+AxdzYMG9+FzwH5JMoc1DmrWW1XXV9VlwG/GUeCAmdR6blXd2c9eADxujmscNJN6fzYwuxUwrqt5ZvK6Bfg/dB/YfjGXxc2UYTF6j66qmwD6n4um65xkb7pPH9fNQW2TrVOtY/BYYNXA/A1929A+VfUr4Dbgd+akuvubSb3zxbrW+qfAl0da0fRmVG+S1yS5ju5N+LVzVNtkzVqT7AEsqao5H9qdqc3GXcBCkOQs4DFDFr1lHbezHfBJ4IiqGsknzdmqdUyGHSFM/rQ4kz5zZT7V0jLjWpP8CbAn8JyRVjS9GdVbVR8EPpjkZcBbgSNGXdgQ09aaZBPgOODIuSpofRgWs6Cq9p9qWZKbk2xXVTf1YbBmin4PBb4EvLWqLhhRqbNS6xjdACwZmH8csHqKPjck2Qx4GHDr3JR3PzOpd76YUa1J9qf7YPGcqrp7jmobZl0f288AJ4y0oqm1at0G2BU4rx8xfQxwepIXVtV356zKBoehRu907vs0cwTwxckdkmwOnAacXFWfncPaJmvWOmYXATsleUL/mB1GV/OgwX04BDin+rOHYzCTeueLZq39UMmHgRdW1bg/SMyk3p0GZv8IuGYO6xs0ba1VdVtVbVtVO1TVDnTng+ZVUABeDTXqiW68/Gy6F+rZwCP79j2Bf+hv/wlwD7BiYFo6H2vt578OrAXuovvU9AdzWOMLgKvpzum8pW97B90vF8BDgM8C1wLfAZ445ue/Ve9e/WP4c+DHwBXzuNazgJsHXqOnz/PH9n3AFX2t5wK7zNdaJ/U9j3l4NZRf9yFJanIYSpLUZFhIkpoMC0lSk2EhSWoyLCRJTYaFJKnJsJAkNf0b/yOuV7XaMNIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from statlearning import plot_coefficients\n",
    "\n",
    "plot_coefficients(stack.meta_regr_, labels = ['OLS', 'Lasso', 'Ridge', 'XGBoost'])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model Evaluation\n",
    "\n",
    "\n",
    "### Original prices\n",
    "\n",
    "To make predictions in the original prices, we exponentiate the y values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Test RMSE</th>\n",
       "      <th>Test R2</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>OLS</th>\n",
       "      <td>14875.800</td>\n",
       "      <td>0.950</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Lasso</th>\n",
       "      <td>14791.127</td>\n",
       "      <td>0.951</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Ridge</th>\n",
       "      <td>14704.635</td>\n",
       "      <td>0.951</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Elastic Net</th>\n",
       "      <td>14791.233</td>\n",
       "      <td>0.951</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Tree</th>\n",
       "      <td>30663.106</td>\n",
       "      <td>0.789</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Bagged Trees</th>\n",
       "      <td>24435.638</td>\n",
       "      <td>0.866</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Random Forest</th>\n",
       "      <td>23904.440</td>\n",
       "      <td>0.872</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SKLearn Boost</th>\n",
       "      <td>27177.019</td>\n",
       "      <td>0.834</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>XGBoost</th>\n",
       "      <td>15955.691</td>\n",
       "      <td>0.943</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Additive Boost</th>\n",
       "      <td>13799.568</td>\n",
       "      <td>0.957</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Stack</th>\n",
       "      <td>13083.293</td>\n",
       "      <td>0.962</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                Test RMSE  Test R2\n",
       "OLS             14875.800    0.950\n",
       "Lasso           14791.127    0.951\n",
       "Ridge           14704.635    0.951\n",
       "Elastic Net     14791.233    0.951\n",
       "Tree            30663.106    0.789\n",
       "Bagged Trees    24435.638    0.866\n",
       "Random Forest   23904.440    0.872\n",
       "SKLearn Boost   27177.019    0.834\n",
       "XGBoost         15955.691    0.943\n",
       "Additive Boost  13799.568    0.957\n",
       "Stack           13083.293    0.962"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "columns=['Test RMSE', 'Test R2']\n",
    "rows=['OLS', 'Lasso', 'Ridge', 'Elastic Net', 'Tree', 'Bagged Trees', 'Random Forest', 'SKLearn Boost', 'XGBoost', 'Additive Boost', 'Stack']\n",
    "results=pd.DataFrame(0.0, columns=columns, index=rows) \n",
    "\n",
    "methods=[ols, lasso, ridge, enet, tree, bag, rf, gb, xbst, abst, stack]\n",
    "\n",
    "for i, method in enumerate(methods):\n",
    "    \n",
    "    if method != stack:\n",
    "        y_pred=np.exp(method.predict(X_test))   \n",
    "        if method == abst:\n",
    "            y_pred=np.exp(lasso.predict(X_test)+method.predict(X_test)) # combining predictions           \n",
    "    else:\n",
    "        y_pred=np.exp(method.predict(X_test.values))\n",
    "        \n",
    "    results.iloc[i,0] = np.sqrt(mean_squared_error(np.exp(y_test), y_pred))\n",
    "    results.iloc[i,1] = r2_score(np.exp(y_test), y_pred)\n",
    "\n",
    "results.round(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Making a Submission on Kaggle\n",
    "\n",
    "Methods from this tutorial would allow one to get a competitive score in the corresponding [Kaggle competition](https://www.kaggle.com/c/house-prices-advanced-regression-techniques). Note that the Kaggle competition is based on predicting the log prices. \n",
    "\n",
    "The next cell shows you how to generate a submission file (note that Kaggle requires inclusion of the Id column, which does not exist in our version of the dataset). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "submission = pd.DataFrame(np.c_[test.index, y_pred], columns=['Id', response])\n",
    "submission['Id'] = submission['Id'].astype(int)\n",
    "submission.to_csv('kaggle_submission.csv',  index=False)"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
